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1 . 0   INTRODUCTION 

As  MOSFET  devices  have  been  shrunk  to  the  submicron 
range  and  a  chip  contains  100,000  transistors  in  the  very 
large  scale  integration  (VLSI)  technology,  trial-and-error 
methodology  to  optimize  such  a  complex  processes  is 
impractical.  Hence  computer  simulation  is  a  cost  effective 
and  time  efficient  alternative,  not  only  providing  an 
appropriate  approximation  before  fabrication,  but  also 
serving  as  a  powerful  tool  in  the  development  of  advanced 
VLSI  technology. 

In  the  design  of  VLSI  processes,  a  series  of 
sophisticated  computer  simulation  programs  from  the  initial 
process  level  to  the  device  level  and  finally  the  circuit 
performance  are  required  [3].  In  the  last  decade,  various 
process,  device,  and  circuit  simulators  have  been  developed 
independently. 

The  reasons  for  integration  are  based  on  two  facts; 
firstly,  Application  Specific  Integrated  Circuits  (ASIC)  have 
emerged  as  a  very  important  IC  product  in  the  market[4]. 
Secondly,  process  and  device  tools  have  been  utilized  to  the 
extent  of  circuit  simulation  since  the  device  features  are 
developed  to  cope  with  the  demand  for  more  highly  integrated 
circuit  like  the  4-Mbit  dynamic  RAM.  Although  the 
benefits  of  simulation  can  be  obtained,  the 


realization  of  using  simulation  programs  is  not  easy  in 
practice.  Three  major  attributes  have  been  identified  to  use 
the  simulation  tools  effectively.  Firstly,  a  good 
understanding  of  the  process,  device,  and  circuit 
characteristics  and  models  being  simulated  is  very  important. 
Secondly,  in  order  to  deal  with  these  sophisticated  and  large 
computer  program,  one  must  be  familiar  with  the  computing 
environment  which  includes  operation  system,  file  management, 
programming  language  in  C  or  FORTRAN  so  that  the  simulation 
programs  can  be  utilized  efficiently.  Finally,  an 
understanding  of  the  numerical  method  is  needed,  either  to 
code  the  mathematical  equations  for  the  modified  models  which 
improve  the  chance  of  obtaining  accurate  results  or  to 
interpret  the  results  when  convergence  problems  arise.  On  the 
basis  of  the  understanding  of  the  above,  all  problems  and 
pitfalls  encountered  while  working  with  numerical  simulation 
programs  can  be  traced  back  and  solved. 

In  this  report,  an  integrated  simulation  for  VLSI  design 
environment  is  proposed.  The  idea  is  to  link  process 
simulator,  device  simulator,  and  circuit  simulator  into  a 
comprehensive  simulation  tool  that  will  improve  the  design 
and  simulation  efficiency  by  automatic  generation  of 
simulation  file  and  automatic  determination  of  device  and 
circuit  characteristics. 

The  architecture  of  this  integration  is  shown  in  Fig   1. 
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FIG.    1.   ARCHITETURE  OF  INTEGRATED  SIMULATORS 


In  the  first  stage,  fabrication  process  is  converted  to 
an  input  file  of  SUPREM.  The  next  stage  is  the  device 
simulator.  Device  geometry  with  the  doping  profile  generated 
by  SUPREM  is  used  to  simulate  the  two  dimensional  device 
characteristics.  The  output  of  PISCES  is  carrier 
distribution,  impurity  distribution,  internal  electric  field, 
electrostatic  potential,  mobility,  and  internal  current 
distribution.  Then,  the  interface  program  is  used  to  generate 
parameters  for  SPICE.  Finally,  based  on  the  parameters 
generated  from  interface  program  and  specified  by  user, 
circuit  performance  is  simulated.  Results  of  the  SPICE, 
PISCES,  and  SUPREM  are  stored  in  the  common  data  base  where 
they  are  employed  by  post  processor  to  print  and  plot  the 
data. 

In   this  report,  most  of  the  capitalized   variables   are 
the  parameters  of  SPICE  which  are  listed  in  appendix  B. 


2.0  PROCESS  SIMULATION 

The  primary  purpose  of  process  simulation  is  to  obtain 
impurity  profiles  in  silicon  which  can  be  used  in  device 
simulation  program  to  theoretically  predict  the  current- 
voltage  characteristics  of  the  fabricated  devices. 

Both  one-  and  two-dimensional  process  simulation 
programs  have  been  developed.  While  the  one-dimensional 
programs  are  well  developed,  two-dimensional  effects  are 
becoming  increasingly  important  as  device  structures  are 
scaled  down  and  lateral  interaction  between  device  regions 
becomes  more  critical. 

2.1  SUPREM:  1-D  Process  Simulator 

SUPREM  stands  for  Stanford  University  PRocess  Engineering 
Modeling  program.  The  program  input  for  SUPREM  is  a 
description  of  the  processing  schedule,  specifying  a  sequence 
of  times,  temperatures,  ambient  and  other  parameters  for 
diffusion,  oxidation,  ion  implantation,  deposition,  and 
etching.  The  output  is  a  one-dimensional  impurity 
concentration  profile  in  the  vertical  direction  for  the 
regions  of  the  silicon  and  some  other  layers,  such  as  silicon 
dioxide  and  polysilicon.  The  block  flow  diagram  for  the 
SUPREM  program,  shown  in  Fig.  2,  indicates  that  the  process 
steps   can  be  simulated  either  individually  or   sequentially, 
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FIG.  2.  FLOW  FOR  THE  SUPREM 


with  the  dopant  profile  predicted  after  one  operation  and 
used  as  the  input  for  the  next  operation.  The  input  file 
resembles  an  actual  process  runsheet  and  consists  of  free 
format  statements.  TITLE,  GRID,  and  SUBSTRATE  cards  are 
always  the  first  three  cards  that  must  appear  at  the 
beginning  for  array  and  parameter  initialization.  After 
impurities  are  introduced,  ion  implantation,  oxidation, 
etching,  deposition,  impurity  predeposition,  and  epitaxial 
growth  are  treated.  The  output  of  the  program  consists  of  the 
one-dimensional  profiles  of  all  the  dopants  present  in  the 
silicon  and  silicon  dioxide  material. 

2.2  Process  Models 

Although  the  fabrication  step  simulation  is  based  on 
several  process  models  in  SUPREM,  only  three  models  for 
technology  steps  which  have  been  used  most  often  and  bring 
out  accurate  results  are  presented  in  this  section. 

2.2.1  Ion  Implantation 

Ion  implantation  provides  a  highly  controlled  technology 
for  introduction  of  the  dopants  into  semiconductor.  The 
simplest  description  of  an  ion  implantation  in  silicon 
dioxide  is  a  symmetric  Gaussian  curve  with  the  two  variables, 
the  projected  range  Rp  and  standard  deviation  RHO.  However, 
the   actual   distributions   of   many   ions,   such   as   boron. 


arsenic,  and  phosphorus,  are  found  to  be  asymmetrical.  Higher 
order  moments,  therefore,  must  be  used  to  construct  range 
distributions  which  are  better  fitted  by  two  half-Gaussian 
profiles,  each  with  different  deviation,  RH01  and  RH02, 
joined  together  at  a  modal  range  Rm.  For  the  joint  half- 
Gaussian  distribution,  the  two  sides  are  given  by 

(2. 1)  C1(y)=Cp*exp[-(y-Rm)2/(2*RH012)] 

(2.2)  C2(y)=Cp*exp[-(y-Rm)2/(2*RH022)] 

where  Cp  is  the  peak  concentration.  SUPREM  can  obtain  the 
three  parameters  from  table  lookup  and  interpolation.  An 
implanted  boron  profiles  in  amorphous  silicon  before 
annealing  has  been  shown  by  W.  K.  Hofker  who  has 
experimentally  determined  the  fourth  moment  to  modify 
Pearson-type  IV  distribution.  Therefore,  an  exponentially 
long  tail  due  to  channeling  can  be  modeled.  The  above 
implantation  model  is  valid  only  for  pre-anneal  data,  since 
enhanced  diffusion  of  implanted  species  occurs  during  any, 
including  a  rapid,  thermal  annealing.  The  details  of  how  to 
model  thermally  annealed  boron  profiles  in  silicon  are 
described  in  Reference  [5],  Fig.  3  shows  different 
characteristics  of  arsenic,  phosphorus  and  boron  impurities 
implanted  with  the  same  energy  and  dose.  Since  boron  is  the 
lightest  atom  among  them,  the  implantation  peak  is  located  at 
the  deepest  position. 
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Fig.  3.  Profiles  of  B,  As  and  P  with  same  energy  and 
dose[ll] . 


2.2.2  Thermal  Oxidation 

The  thermal  oxidation  process  by  which  silicon  dioxide  is 
formed  on  silicon  constitutes  a  major  step  in  the  fabrication 
of  silicon  devices.  It  is  a  direct  reaction  between  atoms  of 
silicon  near  the  surface  of  the  wafer  and  oxygen  supplied  in 
a  high-temperature  furnace.  The  high  quality  of  the  interface 
between  silicon  and  thermally  grown  silicon  dioxide  is 
fundamental  to  the  successful  production  of  MOSFET 
technology.  Three  oxidant  fluxes  involved  in  the  oxidation 
process  are 

(2.3)  Fl=h*(Cl-C0) 

(2.4)  F2=Deff*(CO-C2)/XO) 

(2.5)  F3=k*C2 

where  Fl  is  the  transport  flux  from  gas  ambient  to  the  gas- 
oxide  surface;  h  is  gas  transport  coefficient.  The  transport 
flux  Fl  is  proportional  to  the  difference  between  equilibrium 
concentration  of  oxygen  in  the  oxide,  CO,  and  concentration 
at  the  oxide  surface,  CI.  F2  is  a  diffusion  flux  inside  the 
oxide  layer;  Deff  is  the  effective  diffusion  coefficient;  C2 
is  the  concentration  at  the  interface  between  silicon  and 
oxide;  XO  is  the  thickness  of  the  oxide  layer  already 
present.  F3  is  a  reaction  flux  at  the  silicon/oxide 
interface,   and   k  is  the  surface   reaction   coefficient   for 
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oxidation. 

The  linear-parabolic  oxide-growth  model  in  SUPREM  assumes 
steady-state  oxidant  diffusion  which  implies  that  the  three 
fluxes  are  equal. 

(2.6)  Fl  =  F2  =  F3  =  F 

The  oxide  growth  rate  is  directly  proportional  to  the 
flux  as 

(2.7)  dX0/dt=f/Nl=(k*C2/Nl)/(l+k/h+k*X0/Deff) 

where  Nl  is  the  number  of  oxidant  molecules  incorporated  in  a 
unit  volume  per  unit  time  of  the  oxide  layer.  When 
integrated,  Eq(2.7)  lead  to  the  well-known  linear-parabolic 
growth  relationship  but  only  if  an  initial  oxide  thickness  XI 
is  specified  prior  to  the  oxidation  step  under  consideration 

(2.8)  (X02-Xl2)/B+(XO-Xl)/(B/A)=t 

where  B  is  called  the  parabolic  rate  constant  and  B/A  the 
linear  rate  constant.  Under  relatively  low  dopant 
concentration  conditions,  B  and  B/A  depend  only  on  silicon 
crystal  orientation,  oxiding  ambient  and  temperature. 
Behaviors  of  the  two  rate  are  stored  in  SUPREM  as  constants 
with  default  values  and  are  shown  in  Fig  4. 

Using  high  surface  impurity  concentration,  more  silicon 
vacancies    are   created,   and   the   rate   constants    are, 
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Fig.  4.  (a)Linoar  rate  coefficient  B/A  (b)Parabolic 
rate  coefficient  B[l'^]. 


12 


consequently,  enhanced.  Since  the  impurity  concentration  at 
the  silicon-oxide  interface  changes  due  to  diffusion  and 
segregation  during  oxidation,  enhanced  values  are  calculated 
at  each  time  step  in  SUPREM.  It  is  known  that  the  growth  rate 
is  enhanced  as  much  as  a  factor  of  10  for  oxide  thinner  than 
20  nm  in  a  dry  oxygen  condition.  This  fast  oxidation 
phenomenon  has  attracted  great  deal  of  attention  recently 
because  the  gate  oxide  in  the  present  VLSI  process  is  grown 
within  the  thin  oxide  regime.  SUPREM  III  has  an  empirical 
factor  incorporated  in  the  linear-parabolic  model 

(2.9)  dX0/dt=B/ ( XO+A) +K*EXP< -X0/L ) 

where  the  decay  length  L  is  approximately  independent  of 
temperature  and  K  is  a  singly  activated  function  of 
temperature  with  an  activation  energy  of  2.35  eV  for  < 1 1 1 > 
orientation  and  1.8  eV  for  <100>  silicon. 

2.2.3  Impurity  Redistribution 

Impurity  diffusion  in  silicon  during  high-temperature 
processing  steps  can  be  described  by  the  complete  continuity 
equation 

(2. 10)  dC/dt=d/dx(D*dC/dx)+(-q)/(kT)*d/dx(D*Ci*dPB/dx) 

where  D  is  the  diffusivity,  and  C  and  Ci  are  the  total  and 
electrically   charged  impurity  concentrations,   respectively; 
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PB  is  the  bulk  junction  potential  which  is  expressed  as 
below: 

(2.11)  PB=(k*T/q)*ln(n/ni) 

where  n  and  ni  are  the  extrinsic  and  intrinsic  carrier 
concentrations,  respectively.  The  first  term  in  the 
continuity  equation  represents  classical  concentration 
gradient-driven  diffusion,  including  non-constant 
diffusivity.  The  second  term  incorporates  the  electrostatic 
field-driven  flux. 

SUPREM  utilizes  models  based  on  vacancy  diffusion 
mechanisms  under  non-oxidizing  conditions.  The  intrinsic 
diffusivity  of  an  ionized  impurity  species  is  the  sum  of  the 
dif fusivities  resulting  from  neutral  vacancies  with  an 
opposite  charge.  There  are  four  charge  states  for  vacancies: 
double  negative(=),  single  negative  (-),  neutral  (x),  and 
positive  (+).  Thus,  the  effective  diffusivity  under  non- 
oxidizing  conditions  is 

(2.12)  DN=Dix+  Di"[V"]+Di=[V=]+Di+[V+] 

where  [V~],  [V-],  [V+]  refer  to  the  atom  fractions  of  these 
species  under  extrinsic  conditions  divide  those  under 
intrinsic  condition.  A  boron  atom  as  an  acceptor  is 
negatively  charged  in  the  silicon  lattice  and  it  diffuses 
primarily   with   V+   and   Vx  vacancies.   Thus    the   boron 
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diffusivity  is 

(2.13)  DN(B)  =  Dix+  Di+(p/ni) 

The  diffusion  coefficients  for  boron  and  other  impurities  are 
given  as  default  values  in  SUPREM.  Arsenic  as  a  donor  appears 
to  diffuse  with  Vx  and  V~  vacancies  and  its  diffusivity  is 

(2.14)  DN(As)  =  Dix+  Di~(n/ni) 

The  bandgap  narrowing  effect  due  to  lattice-misfit  strain 
is  also  introduced  in  SUPREM,  which  causes  the  diffusivity  to 

decrease  in  the  heavily  doped  region  that  density  is   greater 

on    _q 
than   5*10   cm   .   Fig.  5  shows  redistribution  profiles   of 

three    implanted  dopants,   given   in   Fig.   3,   where   the 

impurities   were  diffused  separately  at  the  same  temperature 
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Fig.  5.  Redistributed  profiles  of  As,  P,  and  B[ll]. 
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3.0  DEVICE  SIMULATION 

3.1  PISCES-II:  2-D  2-carrier  Poisson  Solver 

The  increased  complexity  of  the  device  structure  is 
necessary  for  optimization  of  the  device  performance,  such  as 
minimizing  the  drain-induced  barrier-lowering  effects, 
enhancing  the  device  reliability,  or  reducing  the  electric 
field  at  the  drain  of  the  MOSFET.  Therefore,  in  the 
development  of  VLSI  MOS  technology,  it  is  essential  to  be 
able  to  simulate  the  electrical  characteristics  of  devices 
which  have  complicated  structures.  The  PISCES  program 
provides  this  capability. 

PISCES-II  is  a  two-dimensional,  two-carrier  semiconductor 
device  modeling  program  which  simulates  the  electrical 
behavior  of  devices  under  either  steady-state  or  transient 
conditions.  It  accepts  data  from  the  SUPREM,  which  provides 
high  accuracy  in  the  impurity  definition,  essential  for 
submicron  device  simulations. 

The  PISCES  program  can  simulate  device  structures  such 
as  MOSFET,  JFET,  MESFET, and  other  non-planar  insulator- 
semiconductor  structures  such  as  the  trench  isolation 
structure.  The  program  can  extract  device  parameters  such  as 
the  threshold  voltage,  subthreshold  slope,  punchthrough 
voltage,  body  effects,  and  electrode  capacitance.  The  program 
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supports  non-uniform  triangular  grids  which  can  be  specified 
through  input  card  or  by  using  an  interactive  grid  generation 
pre-processor  (IGGI).  Further,  PISCES  itself  can  refine  the 
specified  grid  during  the  solution  process.  A  menu-driven 
post-processor  is  available  as  well  for  plotting  solution 
data. 

3.1.2  Basic  Physical  Description  and  Numerical  Algorithm 

To   cover  the  whole  range  of  the  field-effect   transistor 
operation,   PISCES   solves  the  Poisson  equation   and   carriers 
continuity  equations  for  the  electric  potential  and   carriers 
concentrations.  The  Poisson  equation  is 

(3.1)  VPSI=- ( q/EPSILON ) ( p-n+ND+-NA" ) -RHOF/EPSILON 
and  the  continuity  equations  for  electrons  and  holes 

(3.2)  in/   3t=(  VJ)/q+Gn-Rn 

(3.3)  ^p/  j,t=(VJ)/q+Gn-Rn 

where  PSI  is  the  electrostatic  potential  and  EPSILON  is  the 
semiconductor  dielectric  permittivity;  ND+  and  NA~  are 
ionized  donor  and  acceptor  concentrations,  respectively;  RHOF 
is  a  fixed  space  charge  density  which  may  be  present  in 
insulating  material;  n  and  p  are  electron  and  hole 
concentrations,  respectively.  Though  the  Fermi-Dirac 
statistics  should  be  used  in  the  semiconductor,  the  Maxwell- 
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Boltzmann  statistics  is  employed  instead,  because  it  is 
simple  and  still  a  good  approximation.  Then,  electron  and 
hole  density  can  be  expressed  by 

(3.4)  n  =NC*EXP [ ( Ef n-Ec ) /K] =ni *EXP [ q* ( PSI -PH In ) / ( k*T ) ] 

(3.5)  p=NV*EXP[(Ev-Efp)/K]=ni*EXP[q*<PHIp-PSI)/(k*T)] 

where  ni  is  the  intrinsic  carrier  density;  T  is  the  absolute 
temperature,  and  k  is  the  Boltzmann  constant;  PHIn  and  PHIp 
are  the  quasi-Fermi  potentials  of  electrons  and  holes, 
respectively. 

From  Boltzmann  transport  theory,  Jn  and  Jp  can  be 
written  as  function  of  electrostatic  potential  and  the  quasi- 
Fermi  levels  for  electrons  and  holes 

(3.6)  Jn=-q*Un*n*  VPHIn 

(3.7)  Jp=-q*Up*P*7PHIp 

Alternatively,  Jn  and  Jp  can  be  written  as  function   of 
PSI,  n,  and  p,  consisting  of  drift  and  diffusion  components: 

(3.8)  Jn=q*Un*En*n+q*Dn*^"n 

(3.9)  Jp=q*Up*Ep*p+q*Dp*V P 

Where  Un  and  Up  are  the  electron  and  hole  mobilities;  Dn  and 
Dp  are  electron  and  hole  dif fusivities ;  En  =  Ep  =  E  =  -  \7PSI 
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by  neglecting  the  effects  of  band-gap  narrowing  and  assuming 
Boltzmann  carrier  statistics. 

Basically,  the  Poisson  equation  is  solved  numerically 
using  a  five  point  finite-difference  approximation.  The  grid 
structure  used  by  PISCES-II  is  rectangular  in  first  stage  and 
has  nonuniform  spacing  both  in  the  horizontal  and  the 
vertical  directions.  The  grid  spacings  can  be  defined  to  be 
small  in  regions  where  the  potential  changes  rapidly,  such  as 
the  regions  near  the  source,  the  drain,  and  also  near  the 
surface.  A  maximum  of  120  grids  are  allowed  for  the 
horizontal  as  well  as  the  vertical  coordinates.  A  user  can 
specify  total  number  of  grids,  the  minimum  spacing  and 
geometrical  ratio  for  both  x  and  y  directions.  PISCES-II  has 
two  basic  numerical  solution  techniques:  the  Gummel  method 
and  the  full  Newton  method  (sparse  direct  LU  decomposition). 
The  Gummel  method  is  preferred  for  zero/reverse  bias  and  low- 
level  injection  simulation,  while  the  full  Newton  method  is 
good  for  high-level  injection  as  well  as  MOSFET  saturation 
conditions. 

In  Newton's  method,  all  the  variables  PSI,  n,  and  p  are 
solved  for  the  roots  of  each  single  equation  by  Newton- 
Raphson  method  which  is  an  iterative  scheme  to  approximate 
the  root  by  using  successive  tangents.  An  initial  guess  will 
be  given  for  the  unknowns  at  each  node  (PSIo,  no,   po),   and 
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then  the  new  update  value  will  be  calculated  by  solving  the 
linear  system  which  is  formed  by  the  derivative  of  the  matrix 
called  Jacobian  matrix.  This  matrix  can  be  solved  by  Gaussian 
elimination  that  is  by  far  the  most  stable  method.  But 
Newton's  method  takes  time  and  memory  for  two  carriers 
simulations,  for  low  current  solutions  the  gummel  method 
offers  an  attractive  alternative  to  inverting  the  full 
Jacobian. 

The  flow  of  PISCES-II  is  shown  in  Fig. 6.  The  example  of 
pn  diode  and  MOSFET  for  PISCES  will  be  presented  in  Sections 
3.  3  and  3.  4. 
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3.2  Classification  of  Device  Modeling 

A  model  describes  the  behavior  of  a  device  in  a  circuit 
simulation  program.  Ideally,  it  would  be  convenient  to  have 
only  one  model  which  can  fit  all  the  needs.  However, 
different  applications  call  for  different  requirements  of 
device  models.  Basically,  there  are  two  classes  of  device 
models  to  serve  different  applications:  first,  compact 
models,  which  use  the  method  of  Gummel-Poon  or  Ebers-Moll, 
are  based  on  closed  form  solution  to  approximate  device 
equation  for  representing  device  in  circuit  simulation.  The 
corresponding  requirements  for  these  models  are  accuracy, 
simplicity  to  save  computing  time,  understandability,  and 
minimal  number  of  model  parameters  having  direct  visual 
effect  on  device  characteristics.  Second  class  is  the 
theoretical  models  which  use  a  numerical  method,  either 
finite  difference  or  finite  element,  to  solve  nonlinear 
semiconductor  equations  including  Poisson's  equation,  current 
density  equations  for  holes  and  elections,  and  continuity 
equations  for  holes  and  electrons.  This  class  of  models  is 
based  on  the  fundamentals  of  solid-state  physics  for 
predicting  electrical  behavior  of  new  devices.  The 
corresponding  requirements  for  these  models  are  completely 
theoretical  basis  and  understanding  for  the  interpretation  of 
device  physics. 
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Process  technology  of  integrated  circuits  has  been 
improved  to  fabricate  submicron  device  which  encounter  a 
variety  of  physical  effects.  Since  a  theoretical  model 
becomes  difficult  to  manipulate  and  a  compact  model  results 
in  a  loss  of  predictive  capabilities.  A  compromise  is  made  in 
developing  model,  which  is  a  combination  of  physics  based  and 
empirical  equations,  for  circuit  simulation. 

The  device  model  for  circuit  simulation  computes  the 
terminal  currents  and  charges  of  the  device  as  a  function  of 
the  terminal  voltages.  The  terminal  currents  and  charges 
should  be  continuous  functions  of  terminal  voltages  for 
Newton-Raphson  iterations  to  converge.  The  Newton-Raphson 
technique  requires  the  computation  of  the  first  derivatives 
which  are  usually  difficult  to  obtain  analytically.  Hence, 
numerical  derivative  using  either  finite  difference  or  finite 
element  method  is  employed,  but  this  increases  the  model 
computation  time. 

Comparing  to  device  equations,  terminal  capacitance 
equations  are  easy  to  manage,  but  it  is  necessary  to 
formulate  the  charge  equations  if  the  charge  conservation  is 
desired. 

During  model  development,  it  is  important  to  keep  in  mind 
the  modeling  paradox.  Although  more  complex  models  can 
potentially   represent   the  device   characteristics   more 
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accurately,  it  is  more  difficult  to  extract  all  the 
parameters  for  such  complex  models.  Therefore,  if  the  model 
parameters  are  not  specified  properly,  the  device 
characteristics  will  not  be  reproduced  accurately. 
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3.3  Diode  Model  Development 

The  analysis  and  design  of  integrated  circuits  depends 
heavily  on  the  use  of  suitable  models  for  integrated  circuit 
components.  In  order  to  use  models  properly,  the  knowledge  of 
how  the  models  have  been  developed  becomes  important  and 
desirable.  In  this  section,  a  pn  diode  will  be  discussed  by 
using  PISCES. 

In  this  example,  a  basic  simulation  procedure  is 
illustrated.  There  are  usually  three  separate  activities  in  a 
PISCES  simulation.  First,  a  grid  is  created  for  equilibrium 
potential  distributions  which  are  based  on  the  geometry  of 
the  device  structure  and  the  doping.  Second,  solutions  are 
obtained  for  various  bias  points  using  the  grid  which  is 
generated  by  first  step,  with  optional  grid  refinement 
dependent  on  the  solution.  Finally,  solutions  are  analyzed 
and  plotted.  The  three  phases  are  usually  maintained  in 
separate  input  decks.  An  interactive  grid  generation 
preprocessor  (IGGI)  can  be  involved  at  the  first  stage,  and 
the  plotting  post-processor  is  of  use  for  solution  analysis. 

Figure  7  is  an  input  file  for  PISCES  to  generate  and  plot 
the  mesh  for  the  diode  example.  Fig. 8  shows  a  potential 
distribution  for  a  cross  section  of  the  diode.  The  doping 
profile  is  plotted  versus  distance  into  the  device  shown  in 
Fig.   9.   In  Fig.  10  the  input  file  of  the   second   stage   is 
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shown,  in  which  a  new  file  will  be  opened  to  save  the 
calculated  data  into  I-V  or  AC  file.  Fig.  11  shows  a  value  of 
contact  resistor  in  CONTACT  card  and  that  value  is  the  only 
difference  from  the  previous  input  file.  In  Figs.  10  and  11 
the  MODEL  card  specifies  Shockley-Read-Hall  and  Auger 
recombination.  It  also  specifies  the  concentration  dependent 
mobility  from  table  look-up  and  a  lateral  field-dependent 
model  which  has  implemented  for  Si  and  GaAs.  Fig.  12  is  the 
input  file  for  stage  3.  LOAD  card  is  used  to  load  I-V 
information  from  the  log  files  generated  in  the  second  stage. 
The  output  files  have  been  written  in  ASCII  code  rather  than 
the  default  (binary)  by  using  ascii  parameter  in  the  PLOT.  ID 
card.  Fig.  13  shows  I-V  curves  with  the  different  contact 
resistors.  Owing  to  the  lumped  resistance,  actual  contact 
bias  may  differ  from  the  applied  bias  which  can  be  seen  in 
Fig.  14.  In  Fig.  15  the  input  file  for  transient  analysis 
with  step  function  voltage  is  given.  Plots  are  shown  in  Fig. 
16,  17,  and  18.  Transient  analysis  with  ramp  function  voltage 
are  shown  in  Fig.  19  to  Fig.  21.  The  input  file  for  AC 
analysis  is  represented  in  Fig.  22.  The  capacitance  and 
conductance  of  the  pn  diode  are  analyzed  in  Figs.  23  and  24, 
respectively.  Fig.  23  shows  the  diode  capacitive  effects  on 
the  reverse-  and  forward-biased  pn  junction.  There  exists  a 
region  depleted  of  carriers  in  the  reverse  bias  pn   junction 
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in  which  a  parallel  plate  capacitor  is  formed.  When  the 
reverse  voltage  increases,  the  depleted  region  enlarges  and 
the  capacitance  reduces.  When  the  junction  is  forward-biased, 
the  large  number  of  minority  carriers  injected  causes  a  much 
greater  charge  storage  effect,  which  is  known  as  diffusion 
capacitance.  But  when  the  forward-biased  voltage  is  greater 
than  the  built-in  potential  barrier,  the  current  will  flow 
through  the  diode  and  then  the  capacitance  is  too  small  to  be 
measured. 


Title  pn  diode 

$  Generate  and  plot  the  mesh  for  the  diode  example 

$  There  are  two  kinds  of  plot  device  available  in  this 

$  version:  tek4107  and  hp2600  series  plotter. 

option  plotdev=tek4107 

$   The  mesh  card  either  initiates  the  mesh  generation  phase 

$   or  reads  a  previously  generated  mesh. 

mesh    rect  nx=60  ny=3  outf=pn.msh 

$   The  x.m  and  y. m  cards  specify  the  location  of  lines  (1) 

$   of  nodes  (n)  in  a  rectangular  mesh.  The  location  is 

$   interpreted   in  microns.  Since  the  junction   is   at   1.0 

$  micron   specified  by  the  DOPING  card,  the  ratio   0.8   is 

$   selected   to  set  the  space  grid   lines   closely  around 

$   the  junction. 

x.m     n  =  l  1=0  r=l 

x.m     n=20  1  =  1  r=0.8 

x.m     n=60  1=5  r=1.07 

y.m     n=l  1=0  r=l 

y.m     n=3  1=1  r=l 

$   The  region  card  defines  the  location  of  materials  in  a 

$   rectangular  mesh. 

region  num=l  ix. 1=1   ix. h=60  iy. 1=1  iy. h=3  silicon 

$   The  electrode  card  specifies  the  location  of  electrodes  in 

$   in  a  rectangular  mesh. 

elec    num=l  ix. 1=1  ix.h=l  iy. 1=1  iy. h=3 

elec    num=2  ix. 1=60  ix.h=60  iy. 1=1  iy.h=3 

$   The  doping  card  dopes  selected  regions  of  the  device. 

doping   reg=l  n.type  conc=lel5  uniform 

$  JUNC  is  the  location  of  the  junction  (um). 

doping    reg=l   p. type   conc=lel9   x. r=0   y.top=0   y. bot=l 

+         junc=1.0   ratio=l  gauss 

$  Plot. 2d  does  not  work  for  us  so  far. 

plot. 2   bound  no. tic  no. fill  grid  outf =pnl . grid 

$   The  plot. Id  card  plots  a  specific  quantity  along  a  line 

$   segment  through  the  device  (mode  A),  or  plots  an  I-V  curve 

$  of  data  (mode  B) . 

plot.  1   log  dop  abs  a.x=0  b.x=5  b.y=0.5  a.y=0.5  points  ascii 

+       outf=dopl. asc 

$  End  card. 

end 


Fig.  7.  Input  file  of  first  stage  for  diode  device. 
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title  pn  diode 

$  Steady-state  forward  bias  analysis  (w/  contact  resistor) 

option  plotdev=tek4107 

$  Load   the  potential  distribution  generated  in   the   first 

$   stage. 

mesh    in=pn.msh 

$   It  has  a  lumped  resistance  10e+5  ohm/um. 

contac  num=l  resis=le5 

$  SYMB  card  performs  a  symbolic  factorization  in  preperation 

$   for  LU  decomposition  using  Newton  method  in  this  case. 

symb    newton  carr=2 

$  METHOD   provide   parameters   used   to   datermine   the 

$   convergence  of  solution  method  using  automated 

$  Newton-Richardson  in  this  case. 

method   rhsnorm  xnorm  autonr 

models   temp=300  srh  auger  conmob  fldmob  print 

$ 

$   SOLVE  card  instructs  PISCES  to  perform  a  solution  for   one 

$   or  more  specified  bias  points. 

solve    init  ascii  outf =pn2a0-r. asc 

$  Any  I-V  and  AC  data  subsequent  to  the  card  will  be  saved. 

log     outf =IV-r. log 

solve   vstep=0. 1  nsteps=10  elect^l 

solve   vstep=0. 2  nsteps=4  elect=l 

solve   vl=2  ascii  outf =pn2a2-r. asc 

$ 

end 


Fig.  10.  Input  file  of  second  stage  for  diode  device. 
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title  pn  diode 

$   Steady-state  forward  bias  (w/o  resistor) 

mesh  in=pn.msh 

$   If  the   CONTACT   card   is   not   specified,   the    lumped 

$   resistance  is  zero. 

symb  newton  carr=2 

method  rhsnorm  xnorm  autonr 

models  temp=300  srh  auger  conmob  fldmob  print 

$ 

solve  init  ascii  outf =pn2a0-nr . asc 

log  outf=IV-nr. log 

solve  vstep=0. 1  nsteps=10  elect=l 

solve  vstep=0. 2  nsteps=4  elect^l 

solve  vl=2  ascii  outf =pn2a2-nr . asc 

$ 

end 


Fig.  11.  Input  file  of  second  stage  for  diode  without 
contact  resistor. 
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Title  Plot  file 

$  Plot  figures  13  and  14. 

mesh  inf=pn.msh 

$ 

$  FIGURE  13 

Soption   asave=f ig3b. sav 

$  Load  the  I-V  data  generated  in  second  stage. 

plot. Id   inf=IV-nr. log  x. ax=val  y.ax=il  log  abs  points   ascii 

+        outf =pn4. asc 

plot. Id   inf=IV-r.log  x. ax=val  y.ax=il  log  abs  points   ascii 

+        outf =pn4r. asc 

$  FIGURE  14 

plot. Id   inf=IV-r. log   x. ax=val  y. ax=vl  points   min=0   ascii 

+        outf =pn4v. asc 

end 


Fig.  12.  Input  file  of  third  stage  for  plot. 
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title  pn  diode 

$   Transient  analysis  -  voltage  step 

mesh  in^pn.msh 

contac  num=l  resis=le5 

symb  newton  carr=2 

method  rhsnorm  xnorm  autonr 

models  temp=300  srh  auger  conmob  fldmob  print 

$ 

$  Load  the  file  generated  in  the  second  stage. 

load  inf =pn2a2-r. slv 

log  outf=IV-step. log 

solve 

solve  vl=-2  tstep=10e-12  tstop=100e-9 

$  FIGURE  16 

plot. Id  x. ax=time   y.ax-il    x.max=10e-9   points    ascii 

+  outf =  f ig6a. asc 

$  FIGURE  17 

plot. Id  x. ax=time   y. ax=val    x.max=10e-9   points    ascii 

+  outf =fig6bl. asc 

plot. Id  x. ax=time   y. ax=vl    x.max=10e-9    points    ascii 

+  outf =fig6b2 . asc 

$  FIGURE  18 

$  Integration  can  be  performed  in  the  PLOT. ID  card. 

plot. Id  x.  ax=time  y.ax=il  x.max=10e-9  points  integral  ascii 

+  outf=f ig6c. asc 

$ 

end 


Fig.  15.  Input  file  for  transient  analysis. 
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title  pn  diode 

$   Transient  analysis  -  voltage  ramp 

mesh     in=pn.msh 

contac  num=l  resis=le5 

symb     newton  carr=2 

method   rhsnorm  xnorm  autonr 

models   temp=300  srh  auger  conmob  fldmob  print 

$ 

load     inf =pn2a2-r. slv 

log     outf =IV-ramp. log 

solve 

$  Endramp  specifies  the  exact  end  of  the  ramp  in  the  running 

$  time.   It  starts  at  t=tO  and  ends  at  t=endramp. 

solve   vl=-2  tstep=10e-12  tstop=10e-8  endramp=2 . 5e-9 

$ 

option  plotdev=tek4107 

$  FIGURE  20 

plot. Id   x. ax=time   y.ax=il    x. max=10e-9    points    ascii 

+         outf =f ig8b. asc 

$  FUGURE  21 

plot. Id   x. ax=time   y. ax=val    x.max=10e-9    points    ascii 

+         outf=f ig8cl. asc 

plot. Id   x. ax=time   y. ax=vl    x.max=10e-9    points    ascii 

+         outf=fig8c2. asc 

$ 

end 


Fig.   19.  Input  file  for  transient  analysis  with 
ramp  function  voltage. 
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title  pn  diode 
$   AC  analysis 

mesh  in=pn.msh 

symb  newton  carr-2 

method  rhsnorm  xnorm  autonr 

models  temp=300  srh  auger  conmob  fldmob  print 

$  Load  the  file  generated  in  the  second  stagte. 

load  inf =pn2a0-r. slv 

log  acfile=AC. log 

solve  local  vl=-3  ac  freq=le3 

solve  vstep=0. 25  nsteps=12  elect=l  ac  freq=le3 

solve  vstep=0. 1   nsteps=  5  elect=l  ac  freq=le3 

option   plotdev=tek4107 

$  FIGURE  23 

plot. Id  x. ax=vl  y.ax=cl2  points  ascii    outf=f ig9b. asc 

$  FUGURE  24 

plot. Id  x. ax=vl  y. ax=gl2  log  abs  points  ascii  outf=f ig9c. asc 

$ 

end 


Fig.  22.  Input  file  for  AC  analysis 
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3.4  An  Example  of  the  MOSFET  Device  in  PISCES 

A  brief  introduction  of  the  MOSFET  input  file  is 
presented.  Fig.  25  shows  input  deck  of  PISCES  which  is 
automatically  generated  by  MOSMESH  program.  The  structure  and 
doping  information  of  the  MOSFET  for  the  MOSMESH  are  listed 
as  follows: 

lateral  diffusion  constant  =0.8 
poly  gate  width  =  2  um 
spacer  thickness  ;  2  ui 
poly  metal  spacing  =  2  um 
substrate  doping  =  lel5 
gate  oxide  thickness  =  250  ang 
channel  threshold  adjust 

type  doping  -   P 

peak  doping  concentration  =  lel6 

peak  doping  location  =0.1 

dose  =  2el2 
lightly  doped  drain 

type  doping  =  N 

peak  doping  concentration  =  lel8 

peak  doping  location  =  0 

junction  depth  =  0.25  um 
drain  doping 

type  doping  =  N 
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peak  doping  concentration  =  le20 

peak  doping  location  =  0 

junction  depth  =  0. 1  um 
maximum  drain  substrate  reverse  bias  =0.1 
maximum  source  substrate  reverse  bias  =  0 
drain  voltage  =0.1 
starting  gate  voltage  =  -1.0 
sweep  gate  voltage  =  2.0 
gate  bias  increment  =0.1 

After  the  above  values  are  entered  for  the  MOSMESH,  an 
output  file  which  can  be  the  input  file  for  PISCES  is 
produced.  Fig.  26  shows  the  IV  characteristics  of  this 
MOSFET. 
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title  MOSMESH  automatic  mesh  -  Version  8822 


$starting  mesh  card 

mesh  nx=45  ny=22  rect  outf=nmesh 

$x  mesh  locations 

x.mesh  node=l  location=-3 . 1  rat=l 
x.mesh  node=2  location=-3  rat=l 
x.mesh  node=8  location=-l . 59265  rat=1.25 
x.mesh  node=9  location=-l. 37133  rat=1.25 
x.mesh  node=10  location=-l . 15  rat=0.8 
x.mesh  node=ll  location=-l  rat=0. 8 
x.mesh  node=13  location=-0. 72  rat=0. 64 
x.mesh  node=33  location=0. 72  rat=l 
x.mesh  node=35  location=l  rat=1.33 
x.mesh  node=36  location=l . 15  rat=2.51 
x.mesh  node=37  location=l . 37133  rat=1.25 
x.mesh  node=38  location=l . 59265  rat=0. 8 
x.mesh  node=44  location=3  rat=0.8 
x.mesh  node=45  location=3.1  rat=l 

$y  mesh  locations 

y.mesh  node=l  location=-0. 025  rat=l 
y.mesh  node=3  location=0  rat=0.8 
y.mesh  node=4  location=0.  01  rat=l 
y.mesh  node=14  location=0. 35  rat=1.24 
y.mesh  node=19  location=l . 05896  rat=1.25 
y.mesh  node=22  location=2  rat=1.25 

$eliminate  some  nodes  that  aren't  needed 
el  in  y.dir  ix. 1  =  11  ix.h=35  iy. 1=5  iy.h=22 
elim  y.dir  ix. 1=1  ix. h=45  iy. 1=14  iy. h=22 

$region  definitions 

region  num=l  ix. 1  =  1  ix. h=45  iy.  1  =  1  iy.h=3  oxide 

region  num=2  ix. 1=1  ix.h=45  iy. 1=3  iy.h=22  silicon 

$Gate  =  1  Source  =  2  Drain  =  3  Bulk  =  4 
electrode  num=l  ix. 1=11  ix.h=35  iy. 1=1  iy.h=l 
electrode  num=2  ix. 1=1  ix.h=2  iy. 1=3  iy.h=3 
electrode  num=3  ix. 1=44  ix.h=45  iy. 1=3  iy.h=3 
electrode  num=4  ix. 1=1  ix.h=45  iy. 1=22  iy.h=22 

$Doping  Information 

$  junction  depth  is  0. 35um 

doping  reg=2  uniform  x.lef=-3.1  x.rig=3.1 
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+       y.top=0.0  y.bot=2  conc=l . 000000e+15  p. type 

$  Channel  Doping  Region 

doping  reg=2  gauss  cbar=0. 564191  peak=0. 1  cono=2 . 000000e+16 

+       P-type  x.  lef=-3.1  x.rig=3.1  ratio.  la=0.  8 

$  Lightly  Doped  Source 

doping  reg=2  gauss  char=0. 133168   peak=0   conc=l . 000000e+18 

+       n.type  x.lef=-3.1  x.rig=-l  ratio. la=0. 8 

$  Lightly  Doped  Drain 

doping  reg=2  gauss  char=0. 133168   peak-0   conc=l . 000000e+18 

+       n.type  x.  lef  =  1  x.rig=3.1  ratio.la=0.8 

$  Heavily  Doped  Source 

doping  reg=2  gauss  char=0. 0736796  peak=0   conc=l . 000000e+20 

+       n.type  x.  lef=-3.1  x.rig=-1.15  ratio.  la=0.  8 

$  Heavily  Doped  Drain 

doping  reg=2  gauss  char=0. 0736796  peak=0   conc=l . 000000e+20 

+       n.type  x.  lef=1.15  x.rig=3.1  ratio.  la=0.  8 

$  gate  characteristics.  Vd  is  0.100000,  minVg   is   -0  100000 

$   maxVg  is  2.000000 

symb  gummel  carriers^l  electrons 

method  iccg  damped 

models  srfmob 

$  change  the  below  line  if  another  gate  material  is  used 

contac  num=l  n.poly 

$  solve  for  the  initial  step.   save  in  file  initial 

solve  initial  vl=-0. 100000  outf ile=initial 

$  save  I-V  information  in  gatelV. log 

log  ivf ile=gateIV. log 

regrid  potential  step=0. 2 

symb  gummel  carriers=l  electrons 

method  iccg  damped 

solve   vl=-0. 100000  v3=0. 100000  vstep=0. 100000   nsteps=21 

+       electrode=l 

$  plot  statement  in  ascii  format  in  file  iv.ascii 

plot. Id  x. axis=vl  y. axis=i3  ascii  outf ile=iv. ascii 

end 


Fig.  25.  A  input  file  of  M0SFET  for  PISCES. 
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4.0  CIRCUIT  SIMULATION 

4.1  SPICK:  Circuit  Simulator 

SPICE,  which  stands  for  simulation  program  with 
integrated  circuit  emphasis,  is  a  general-purpose  circuit 
simulation  program  developed  at  the  University  of  California 
at  Berkeley  for  nonlinear  dc,  nonlinear  transient,  and  linear 
ac  analysis.  Circuits  can  contain  resistors,  capacitors, 
inductors,  mutual  inductors,  independent  voltage  and  current 
sources,  four  types  of  dependent  sources,  transmission  lines, 
and  the  four  most  common  semiconductor  devices:  diodes,  BJTs, 
JFETs,  MESFET,  and  MOSFET. 

SPICE  has  built-in  models  for  the  semiconductor  devices, 
and  the  user  need  specify  only  the  pertinent  model  parameter 
values.  The  model  for  the  BJT  is  based  on  the  integral  charge 
model  of  Gummel-Poon;  however,  if  the  Gummel-Poon  parameters 
are  not  specified,  the  model  reduces  to  the  simpler  Ebers- 
Moll  model.  In  either  case,  charge-storage  effects,  ohmic 
resistances,  and  a  current-dependent  output  conductance  can 
be  included.  The  diode  model  can  be  used  for  either  junction 
diodes  or  Schottky  barrier  diodes.  The  JFET  model  is  based  on 
the  FET  model  of  Shichman-Hodges .  Four  MOSFET  models  are 
implemented:  M0S1  is  described  by  a  square-law  I-V 
characteristic,   M0S2  is  an  analytical  model,  and  M0S3   is   a 


f.i: 


semiempirical  model.  M0S4  is  the  Berkeley  Short-channel  IGFET 
Model  (BISM).  M0S2 ,  MOS3,  and  M0S4  include  second-order 
effects  such  as  channel-length  modulation,  subthreshold 
conduction,  scattering  limited  velocity  saturation,  small- 
size  effects,  and  charge-controlled  capacitances[ 1 ] . 

SPICE  is  node  voltage-oriented,  so  any  node  voltage  can 
be  requested.  Element  currents  flowing  through  independent 
voltage  sources  can  also  be  requested.  Tabular  lists  and 
printer  plots  are  available. 

4.2  Semiconductor  devices  in  SPICE 

The  models  for  the  five  semiconductor  devices  included 
in  SPICE  require  many  parameter  values.  Moreover,  many 
devices  in  a  circuit  often  are  defined  by  the  same  set  of 
device  model  parameters.  For  these  reasons,  a  set  of  device 
model  parameters  is  defined  on  a  separate  .MODEL  card  and 
assigned  a  unique  model  name.  The  device  element  cards  in 
SPICE  then  reference  the  model  name.  This  scheme  alleviates 
the  need  to  specify  the  same  model  parameters  on  each  device 
element  card. 

Each  device  element  card  contains  the  device  name,  the 
nodes  to  which  the  device  is  connected,  and  the  device  model 
name.  In  addition,  other  optional  parameters  can  be  specified 
for  each  device:  geometric  factors  and  an  initial  condition. 

The  value  of  input  parameters  for  the  five  semiconductor 
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device  models  will  be  calculated  to  determine  the  accuracy  of 
the  models.  These  input  parameters  should  be  related  to  the 
particular  process  used  at  each  manufacturing  site.  In  the 
next  section,  MOSFET  models  will  be  emphasized  and  discussed. 


54 


4.3  MOSFET  Model 

The  metal -oxide-semiconductor  field  effect  transistor 
(MOSFET)  is  the  most  important  device  for  for  VLSI  circuits. 
Figure  27  shows  the  structure  of  a  MOSFET  in  which  the  region 
covered  by  the  gate  and  thin  oxide  is  the  channel,  and  its 
length  is  L.  This  length  is  indicated  in  the  SPICE  input  deck 
of  the  MOS  transistor.  Leff  is  the  effective  channel  length 
which  is  defined  as  the  distance  between  two  diffused 
regions.  The  width  of  the  channel,  H,  is  the  width  of  area 
covered  by  the  thin  oxide  between  the  thick  field  oxide. 

4.3.1  The  MOS  Capacitor 

The  MOS  capacitor  is  not  only  the  basis  of  the  MOS 
transistor,  but  is  a  useful  nonlinear  device  whose 
nonlinearity  is  fundamental  to  a  number  of  circuit 
techniques.  The  MOS  capacitor  combines  the  characteristics  of 
the  parallel  plate  capacitor,  Cox,  and  the  pn  junction  diode. 
Since  the  charges  on  a  plate  will  be  mirrored  by  charges  on 
another  plate,  Cox  is  independent  of  the  applied  voltage.  The 
pn  junction  capacitor,  Csi,  which  has  been  discussed 
previously  is  voltage  dependent.  The  silicon  surface 
represents  an  equipotential  plane  that  acts  virtually  as  a 
common  plate  for  both  capacitors  in  series.  In  order  to 
describe  the   latter  capacitance  by  an   equivalent  parallel 
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plate  capacitor,  we  must  determine  the  effective  depth  below 
the  silicon  surface  at  which  an  imaginary  bulk  plate  should 
be  placed.  This  effective  depth  is  referred  to  as  the  Debye 
Length. 

Since  we  do  not  contact  the  silicon  surface,  we  induce  a 
voltage  variation  between  the  silicon  interface  and  the 
silicon  bulk  by  varying  the  voltage  between  the  gate  and  the 
silicon  bulk.  When  a  voltage  variation  is  positive  with 
respect  to  p-type  substrate,  a  positive  charge  is  induced  at 
the  metal-oxide  interface.  The  positive  charge  at  the  metal- 
oxide  repels  holes  from  silicon  surface,  leaving  behind 
negatively  charged,  unneutral ized  and  immobile  acceptor  ions. 
Note  that  the  charge  of  acceptor  ions  balances  the  positive 
charge  induced  at  the  metal-oxide  interface  and  the  majority 
carriers,  holes,  follow  the  signal  changing  at  the  terminal. 
The  difference  in  the  depth  of  distribution  of  free  carriers 
between  the  metal  and  silicon  lies  in  the  concentration  of 
the  free  carriers.  The  concentration  of  electrons  in  metal  is 
in  the  order  of  10  cm  and  in  silicon  is  roughly  in  the 
order  of  10  cm  .  Consequently,  the  depth  of  distribution 
of  free  carriers  that  follow  the  signal  is  considerably 
smaller  in  metal  than  in  silicon.  The  number  of  free  carriers 
that  follow  the  voltage  decreases  exponentially  with 
increasing  depth.  Therefore,  deep  in  the  silicon  bulk  the 
free   carriers  which  are  shielded  by  the  free  carriers   lying 
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above  them  do  not  "feel"  the  voltage  variation.  The  depth 
below  the  surface  at  which  the  number  of  free  carriers  that 
follow  the  voltage  signal  has  decreased  to  37%  of  the  free 
carrier  concentration  is  the  Debye  Length.  This  results  in  a 
capacitance  per  unit  area  that  is  less  than  Cox,  and  is  the 
series  connection  of  these  two  capacitors.  Therefore, 

(4.1)     1  /  Cmos  =  1  /  Cox  +  1  /  Csi 

The  capacitance  of  a  MOS  capacitor  in  inversion  region  is 
much  more  complicated  than  it  was  in  depletion  and 
accumulation.  Because  the  electrons  which  are  observed  in 
equilibrium  at  the  inversion  layer  come  from  the  thermal 
generation  of  electrons  at  the  oxide-silicon  interface.  This 
process  is  slow  which  is  in  the  order  of  millisecond.  Thus  if 
we  measure  the  capacitance  of  the  MOS  structure  in  inversion 
very  slowly,  only  the  capacitance  value  of  Cox  is  observed. 
This  is  because  the  thermal  generation  of  carriers  will  short 
out  the  depletion  layer  capacitance.  But  when  the  applied  ac 
signal  is  changed  rapidly  while  the  dc  bias  voltage  is  varied 
slowly,  the  inversion  layer  cannot  respond  to  the  ac  applied 
signal.  The  number  of  charges  in  the  silicon  space-charge 
layer  is  modulated  instead  of  the  movement  of  the  holes  at 
the  far  edge  of  the  depletion  region.  The  capacitance  then 
corresponds   to   the   series   combination   of   the   oxide 
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capacitance  and  the  depletion  region  capacitance. 

Since  the  depletion  region  reaches  a  maximum  width  when 
the  device  goes  to  strong  inversion,  the  depletion  region 
capacitance  remains  constant  as  the  bias  voltage  is  increased 
further.  This  high  frequency  C-V  curve  is  shown  by  the  solid 
line  in  Fig  28. 

4.3.2  MOS  Transistor  Model 

In  SPICE3A7  there  are  four  different  MOS  models 
available  to  the  users.  The  Level  1  model  is  the  simple 
Shichman-Hodge  model  which  has  been  found  necessary  for 
checking  out  the  correctness  of  hand  calculations  and  for 
understanding  or  developing  new  circuit.  The  Level  2  model  is 
an  analytical  one-dimensional  model  which  incorporates  most 
of  second  order  effects  of  small  size  devices.  The  Level  3 
model  is  a  semi-empirical  model  described  by  a  set  of 
parameters  which  are  defined  by  curve-fitting  rather  than 
physical  background.  The  Level  4  model  is  the  Berkeley  Short- 
channel  IGFET  Model  (BSIM)  which  is  used  in  the  design  with  a 
process  characterization  system  that  provides  all  the 
parameters.  There  are  no  default  values  for  the  parameters 
and  leaving  one  parameter  out  is  considered  an  error.  Thus, 
it  is  necessary  for  the  circuit  designer  to  know  what  are  the 
equations  governing  the  behavior  of  the  MOSFET. 
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4.3.3  Surface  Behavior  of  MOSFKTS 

The  first  step  in  studying  the  theory  of  the  MOSFET  is 
a  behavior  of  the  surface  of  a  semiconductor  under  the 
influence  of  an  electric  field.  This  electric  field,  which  is 
produced  by  the  voltage  applied  between  the  gate  and  the 
substrate,  is  perpendicular  to  the  oxide-semiconductor 
interface. 

The  reference  condition  of  the  surface  is  when  the 
semiconductor  has  the  same  carrier  concentrations  at  the 
surface  as  at  the  substrate.  In  other  words,  the 
concentrations  of  electrons  and  holes  are  not  disturbed  at 
the  silicon  surface  when  no  bias  voltage  is  applied  to 
structure.  This  state  is  called  the  flat-band  condition  and 
is  achieved  by  the  application  of  a  voltage  equal  to  flat- 
band  voltage  between  the  metal  and  the  silicon, 

(4.2)  VFB  =  PHIms  -  QO  /  Cox 

where  QO  is  positive  charge  at  the  oxide-silicon  interface, 
and  therefore  it  requires  a  negative  charge  to  neutralize  its 
effect  and  Cox  is  the  capacitance  per  unit  area  of  the  thin 
oxide  layer.  When  VGB  =  VFB,  the  carrier  concentration  is 
constant  in  the  semiconductor;  in  other  words,  it  is  equal  to 
NA.  The  condition  VGB  <  VFB  increases  the  negative  charge  on 
the  gate  and  increases  the  concentration  of  the  holes  near 
the  surface   in  order  to  balance  the  gate  charge.   Thus   the 
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concentration  of  the  p-type  carriers  at  the  surface  is 
greater  than  at  the  substrate,  and  this  state  is  called 
an  accumulation.  When  VGB  >  VFB.  the  holes  are  pushed  away 
from  the  surface  so  that  the  negative  charge  of  the  fixed 
ions  restores  the  balance  with  the  gate  charge.  The  carriers 
concentration  near  the  surface  is  less  than  that  in  the  bulk, 
and  this  state  is  called  a  depletion.  The  thickness  WO  of  the 
depleted  region  is  shown  in  Fig.  29  and  the  equation  is 
listed  below: 

(4.3)  WO  =  SQRT(2*EPSIL0Nsi*PB/q*NA) 

where  PB  is  the  potential  across  the  depleted  region.  It  is 
interesting  to  know  how  the  carriers  vary  in  the 
accumulation,  depletion,  and  inversion  regions  when  a  voltage 
is  applied.  A  solution  for  this  electrostatic  potential  as  a 
function  of  y  direction,  which  is  perpendicular  to  the 
surface  and  goes  into  the  semiconductor,  allows  a  complete 
description  of  the  carrier  concentration.  The  electrostatic 
potential  is  defined  as 

(4.4)  PSI=[Ei(y)-Ei(bulk)]/q 

The  electrostatic  potential  evaluated  at  the  surface, 
y=0,  will  be  called  the  surface  potential,  PHI.  For  a 
grounded  bulk,  the  potential  and  the  electric  field  are 
related  by 
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(4.5)  E  =  -  V  PSI 

To  obtain  a  solution  for  PHI,  Poisson's  equation  must  be 
solved.  In  one  dimension,  Poisson's  equation  is 

(4.6)  d  PSI/d  y   =  -  RHO/EPSILONs 

where  RHO,  the  total  space-charge  density,  is  defined  as 

(4.7)  RHO  =  q*(Pp  -  Np  +  ND+  -  NA+) 

The  majority  carrier  concentration,  Pp,  and  the  minority 
carrier  concentration,  Np,  can  be  expressed  as 

(4.8)  Pp=ni*exp(q*PSI/kT) 

(4.9)  Np=ni*exp(-q*PSI/kT) 

where  ni  is  the  intrinsic  carrier  concentration  carrier 
concentration.  In  the  bulk,  phi  =0  and  RHO  =  0,  hence 

(4.10)  ND+  -  NA~  =  Npo  -  Ppo 

where  Npo  and  ppo  are  the  equilibrium  densities  of  electrons 
and  holes,  respectively,  in  the  bulk  of  the  semiconductor. 
The  densities  of  the  surface  are 

(4.11)  Ns  =  Npo  *  exp(q*PHI/(k*T) ) 

(4.12)  Ps  =  Ppo  *  exp(q*PHI/(k*T) ) 
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Therefore,  Poisson's  equation  may  be  expressed  as 

(4. 13)    PSI=-q/RHO*[Ppo*( (exp( -z )-l)-Npo*( (exp(z)-l)] 

where  z  is  equal  to  q*PSI/(k*T).  Integrating  Eq.  (4.13)  from 
the  bulk  toward  the  surface,  transforming  it  into  a  form 
relating  potential  PSI  to  the  electric  field,  and  using 
Gauss's  law,  the  space  charge  per  unit  area  in  the 
semiconductor  can  be  expressed  as 

(4.  14)  Qs=-EPSIL0Nsi*E=-1.414*EPSIL0Nsi/(m*LD)*F(PSI) 

where  LD  is  the  extrinsic  Debye  length  for  holes,  and  m  is 
equal  to  q/(k*T). 

Figure  30  shows  the  variation  in  the  space  charge  as   a 
function  of  the  surface  potential.  Notice  the  sharp   increase 
in   charge  as  surface  increases  past  2*PB.  Hence,  the  onset   of 
strong  inversion  is  defined  as  occurring  when 

(4. 15)    PHI=2*PB=2*(k*T/q)*ln(NA/ni) 

4.3.4  Level  1  Model  in  SPICE 

The  level  1  MOSFET  model  used  in  SPICE  is  developed  by 
Shichman  and  Hodges[7].  In  this  model,  the  equations  used  for 
linear  and  saturation  regions  are  listed  below: 

linear  region 

For  VGS  >  VTH  and  VDS  <  VGS  -  VTH: 
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(4. 16)  IDS=KP*[W/(L-LD)]*(VGS-VTH-VDS/2)*VDS*(1+LAMBDA*VDS) 

(4. 17)  VTH=VT0+GAMMA(SQRT(2*PB-VBS)-SQRT(2*PB) 

saturation  region 

For  VGS  >  VTH  AND  VDS  >  VGS  -  VTH : 

(4. 18)  IDS=(KP/2)[W/(L-2*LD)]*(VGS-VTH)"2*(1+LAMBDA*VDS) 

where  W  and  L  are  the  values  of  the  length  and  width  of  the 
channel  specified  in  the  SPICE  input  card.  The  amount  by 
which  the  gate  electrode  overlaps  the  source  and  drain 
regions  must  be  subtracted  from  the  nominal  channel  length. 
Thus, 

(4. 19)  Leff  =  L  -  2  *  LD 

The  term  1  +  LAMBDA  *  VDS  introduced  in  the  model  is  an 
empirical  correction  of  the  conductance  in  the  saturation 
region.  There  are  five  parameters  which  refer  to  the 
electrical  behavior  of  the  MOSFET  in  the  this  model:  KP,  VTO, 
GAMMA,  PHI,  and  LAMBDA. 

They  can  be  specified  directly  in  the  .MODEL  card,  or 
they  can  be  calculated  from  physical  parameters,  using  the 
following  equations: 
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(4.20)  KP  =  u  *  Cox 

(4.21)  GAMMA  =  SQRT( 2*EPSIL0Ns*q*NA) /Cox 

(4.22)  2  *  PB  =  2*(kT/q)*ln(NA/ni) 

Where  Cox  =  EPSILONs  /  Tox.  In  case  of  conflict,  the  value 
of  the  electrical  parameter  is  read  from  the  input  ,  but  they 
will  not  be  computed  from  the  above  equations. 

4.3.5  Equations  for  the  Level  2  Model 

To  obtain  a  better  model  for  IDS,  it  is  necessary  to 
consider  QB  which  will  be  affected  by  the  voltage  in  the 
channel.  According  to  Meyer  model,  we  have 

(4.23)  dQB=W*dx*GAMMA*Cox*SQRT(2*PB-VBS+VC(x)) 
Therefore,  a  new  equation  for  the  current  can  be  obtained  as 

(4.24)  IDS=BETA{(VGS-VFB~2*PB-VDS/2)*VDS-2/3*GAMMA[(VDS-VBS 
+2*PB)"1.5]} 

Where  BETA  is  u*Cox*(W/Lef f ) . 

(4. 25)  VD(sat)=VGS-VFB-2*PHI+GAMMA"2*[l-SQRT(l+2/GAMMA* 

(VGS-VFB)) 

This  model  has  validity  limits  that  can  be  easily  reached 
with  the  present  technology.  The  method  chosen  to  cope  with 
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this  problem  is  to  adapt  the  model  by  applying  semi -empirical 
corrections  to  the  basic  equations  to  develop  new  and  more 
accurate  model  rather  than  start  again  form  the  basis  of  the 
theory. 

4.3.6  Level  2  Model  in  SPICE 

This  section  summarizes  the  basic  equations  of  Meyer's 
model  as  they  have  been  introduced  in  SPICE.  The  corrections 
which  are  made  to  simulate  effects  not  provided  for  in  the 
theory  of  the  basic  model  in  the  previous  section  will  be 
described. 

The  threshold  voltage  can  be  calculated  from  physical 
parameters  as  follows 

(4.26)  VT0=PHIms-q*NSS/Cox+2*PB+GAMMA*SQRT(2*PB) 

Where  PHIms  is  the  potential  difference  between  gate  and 
substrate.  The  current  in  the  linear  region  is  calculated 
through  Eq. (4.24),  where  the  corrective  term  of  the  channel 
length  is  also  present.  Therefore, 

(4.27)  IDS=BETA*{ ( VGS-VFB-2*PB-VDS/2 ) *VDS-2/3*GAMMA[ ( VDS 
-VBS+2*PB)"1.  5]} 

Where  BETA  is  [KP/( 1-GAMMA*VDS) ] *[W/< L-2*LD) .  The  current  in 
the  saturation  region  is 
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(4.28)    IDS=ID(sat)*[l/(l-GAMMA*VDS) 

Where  ID(sat)  is  calculated  from  last  equation  at  VDS  = 
VD(sat),  and  VD(sat)  is  from  Eq.  (4.25).  These  equations  give 
better  result  than  level  1  model,  but  they  are  still  not 
sufficient  for  a  good  agreement,  with  experimental  data,  even 
for  long  channel  simulation.  Therefore,  the  following  will 
show  the  modification  of  the  above  equations  to  improve 
precision  for  long-  and  wide-channel  MOSFET  of  level  2  model. 
Fig.  31  shows  how  well  the  level  1  model  and  level  2  model 
simulate  the  measured  data. 

Effect  of  the  gate  electric  field  on  mobility 

The  constant  mobility  has  been  assumed  with  the  applied 
voltage  so  far.  A  reduction  in  mobility  with  an  increase  in 
the  gate  voltage  has  been  observed.  In  order  to  simulate  this 
effect,  a  variation  of  the  parameter  KP  is  introduced  in 
SPICE: 

(4.29)   KP,=KP*{(EPSILONs/EPSILONox)*[UCRIT*Tox/(VGS-VTH- 
UTRA*VDS)]"UEXP 

Where  UCRIT  is  the  gate-to-channel  critical  field  parameter 
for  mobility  degradation  to  determine  when  the  mobility 
begins  to  decrease  and  UTRA  is  the  transverse  field 
coefficient   chosen  between  0  and  0.5  (dimensionless)  due   to 
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Fig.  31.  (a)  Output  of  level  1  model  in  SPICE  with 

different  KP  to  simulate  the  measured  data, 
(b)  Output  of  level  2  model  in  SPICEC15], 
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drain  voltage.  The  use  of  this  formula  has  been  proven  a  good 
agreement  between  SPICE  and  experimental  data. 

Subthreshold  conduction 

In  the  basic  model  implemented   in   SPICE,   the  drift 

current   is   calculated  only  when  the   surface  potential  is 

equal  to  or  greater  than  2*PB.  In  reality,  a  concentration  of 

electrons   near  the  surface  exists  even  for  VGS   <   VTH,  and 

therefore  there   is   a  current  when  the  surface   is   not  in 

strong  inversion.  This  current  is  due  to  a  diffusion  of 
charge  carriers  between  the  source  and  the  channel. 

(4.30)  IDS=ION*exp[(VDS-VON)*q/(n*k*T)] 

where  ION  is  the  current  in  the  weak  inversion  for  VGS=V0N, 
and  VON,  shown  in  Fig.  32,  is  defined  as 

(4.31)  VON=VTH+z*k*T/q 

where  z=l+(q*NFS+Cd) /Cox.  NFS  is  defined  as  the  number  of 
fast  surface  states  which  also  determines  the  slope  of  IDS 
vs.  VGS  characteristics  in  the  weak  inversion,  and  Cd  is  the 
capacitance  associated  with  the  depleted  region  and  is  given 
by: 

(4.  32)   Cd=dQB/dVDS=GAMMA/[2*SQRT(2*PB-VBS)] 

It  is  clear  that  the  model  introduces  a  discontinuity   in 
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the  derivative  and  therefore  the  simulation  of  the  transition 
region  between  strong  and  weak  inversion  is  not  very 
accurate. 

Variation  of  channel  length  in  the  saturation  region 

An  empirical  equation  using  the  parameter  LAMBDA  has  been 
used   to  calculate  the  conductance  in  the  saturation   region. 
The  level  2-  model  also  offers  a  physical  model  to  compute  the 
channel  length  in  saturation: 

(4.33)  L'=Leff(l-LAMBDA*VDS) 

If  the  parameter  LAMBDA  is  not  specified  in  the  .MODEL 
card,  it  will  be  calculated  using  the  following  equation: 

(4.34)  LAMBDA= ( Lef f -L ' )/(Leff*VDS) 

Where, 

(4.35)  Leff-L'=XD*{(VDS-VD(sat))/4+sqrt[l+(VDS-VD(sat))/4]} 
and 

(4.36)  XD=SQRT[2*EPSIL0Ns/(q*NA)] 

This  model  shows  a  correct  dependence  of  the  conductance 
in  the  saturation  region  on  Leff.  Thus,  if  the  model  with 
constant  LAMBDA  is  used,  the  transconductance  KP  does  not 
vary  with  Leff  when  VDS  >  VD(sat). 
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Effect  of  channel  length  on  threshold  voltage 

Experimental  data  has  shown  that  when  the  channel  length 
is  small  enough  to  compare  with  the  thickness  of  the  depleted 
region,  the  channel  dimensions  W  and  L  will  have  an  effect 
on  the  threshold  voltage.  This  phenomenon  can  be  explained  in 
two  ways  from  the  following  equation 

(4.37)  VTH=VFB+2*PB-QB/Cox 

A  reduction  in  the  charge  QB  due  to  the  source  and 
drain  depleted  region  below  threshold,  or  an  increase  in  the 
channel  surface  potential  due  to  the  effect  of  the  voltage 
VDS.  This  effect  can  be  introduced  by  modifying  the  value  of 
GAMMA  as 

(4.38)  GAMMA'=  GAMMA*{l-[XJ/( 2*Lef f ) ] *[SQRT( 1+2*WS/XJ) 

+SQRT( l+2*WD/XJ)-2)] 

where  HS  and  WD,  shown  in  Fig.  33,  are  the  widths  of  the 
depleted  regions  of  source  and  drain,  respectively,  and  are 
defined  as 

(4.39)  WS=XD*SQRT(2*PHI-VBS) 

(4.40)  WD=XD*SQRT(2*PHI-VBS+VDS) 

Where  XD=SQRT[2*EPSIL0Ns/( q*NA) ] .  In  practice,  a  correct 
analytical  model  of  two-dimensional  effects  on  VTH  is  very 
tedious,   but  it  is  possible  to  approximate  a  good  model   for 
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threshold  voltage  of  a  short-channel  MOSFET  with  two 
dimensional  numerical  model.  This  type  of  analysis  is  time- 
consuming  for  a  circuit  simulation  program,  but  it  can  be 
used  in  the  study  of  physical  behavior  of  short-channel 
MOSFET. 

Effect  of  speed  limit  of  the  carrier 

The   calculation   of   the  saturation   voltage   using  Eq. 

(4.25)   is   based  on  the  hypothesis  that  the   charge   in  the 

channel  is  zero  after  pinch  off  point  or  near  the  drain,  but 

it   actually  has  to  maintain  a  minimum  concentration   of  the 

carriers   to   sustain   the   saturation   current.  This 

concentration  depends  on  the  speed  at  which  the  carriers  are 
moving. 

The   electric  field  between  the  drain  and  the  pinch  off 

point   can  be  sufficiently  high  to  drift  the  carriers  at  the 

speed  limit,  which  is  the  maximum  speed  allowed  for  the 
scattering  effect  within  the  crystal  lattice.  The  VMAX  can  be 
obtained  by 

(4.41)   VMAX=ID(sat)/(W*QI) 

where  QI  is  the  mobile  charge  in  the  channel  for  VDS=VD(sat). 

The  model  discussed  so  far  provides  good  results  in   the 

simulation  of  MOSFET  with  a  minimal  channel  length  of  about  4 

to  5  um,  but  this  model  has  a  discontinuity  in  the  derivative 
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at  the  boundary  between  the  saturation  region  and  linear 
region,  and  also  between  subthreshold  and  linear  region.  This 
phenomenon  results  in  a  less  precise  calculation  of  the 
conductance  and  also  the  cause  of  difficulties  in  the 
convergence  of  the  Newton-Raphson  algorithm. 

4.3.7  Gate  Capacitance 

SPICE  uses  a  gate  capacitance  model  similar  to  that 
proposed  by  J.  E.  Meyer [16].  In  this  simple  model,  the 
charge-storage  effect  is  represented  by  three  nonlinear  two- 
terminal  capacitors:  CGB,  CGS  and  CGD.  The  equations  of  this 
model  are  listed  below  and  shown  in  Fig.  34. 

Accumulation  region 

For  VGS  <  VON  -  2*PB: 

(4.42)  CGB  =  Cox  +  CGBO  *  Leff 

(4. 43)  CGS  =  CGSO  *  W 

(4. 44)  CGD  =  CGDO  *  W 

Depleted  region 

For  VON  -  2  *  PB  <  VGS  <  VON: 

(4.4b)  CGB  =  Cox*  (V0N-VGS)/(2*PB)+CGB0*Leff 

(4.46)  CGS  =  2/3*Cox*[(V0N-VGS)/(2*PB)+l]+CGS0*W 

(4. 47)  CGD  =  CGDO  *  W 

7  b 
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Fig.  34.  Meyer's  model  of  the  capacitance[ 15 ] 
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Saturation  region 

For  Von  <  VGS  <  VON  +  VDS : 

(4. 48)  CGB  =  CGBO  *  Leff 

(4.49)  CGS  =  2/3*Cox  +  CGSO*  W 

(4. 50)  CGD  =  CGDO  *  W 

Linear  region 

For  VGS  >  VON  +  VDS: 

(4.51)  CGB  =  CGBO  *  Leff 

(4.52)  CGB  =  Cox{l-[(VGS-VDS-VON)/(2VGS-2VON-VDS)]"2}+CGSO*W 

(4.53)  CGD  -    Cox{l-[(VGS-VDS-VON)/(2VGS-2VON-VDS)]"2}+CGDO*W 

The  voltage  VON  is  calculated  from  Eq.  (4.31).  If  the 
parameter  NFS  is  specified  in  the  MODEL  card,  then  VON  is 
equal  to  VTH. 

CGBO,  CGSO,  and  CGDO  are  the  overlap  capacitances  among 
the  gate  electrode  and  the  other  terminals  outside  the 
channel  region.  The  model  introduced  in  SPICE  differs  from 
that  developed  by  Meyer  because  of  the  discontinuity  between 
the  capacitance  with  the  accumulation  surface  and  the 
depleted  surfaces  causes  nonconvergence  in  the  Newton-Raphson 
algorithm. 
4.3.8  Junction  Capacitance 

The  capacitance  of  the  diffused  regions  of  the  source  and 
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drain  is  simulated  with  the  pn  junction  model.  The  total 
capacitance  of  a  diffused  region  is  calculated  from  the  sum 
of  an  area  and  a  perimeter  capacitance: 

(4.  54)  CBS=(CJ*AS)/(1-VBS/PB)-MJ+(CJSW*PS)/(1-VBS/PB)-MJSW 

(4.55)  CBD=(CJ*AD)/(1-VBS/PB)~MJ+(CJSW*PD)/(1-VBS/PB)~MJSW 

where  CJ  is  the  bulk  capacitance  at  zero-bias  voltage  per 
square  meter  and  CJSW  is  zero-bias  perimeter  capacitance  per 
meter,  and  PB  is  the  junction  potential  that  can  be  put  in 
the  .MODEL  card  or  calculated  from  the  physical  parameters  as 

(4.53)  PB=Eg/2+(kT/q)*ln(NB/ni) 

The  exponential  factors  have  the  following  default 
values:  MJ  =  0.5,  for  the  step- junction  approximation,  and 
MJSW  =  0.33,  for  the  linearly  graded  approximation.  The 
junction  capacitances  are  very  important  not  only  for  the 
correct  description  of  the  circuit  to  be  simulated  but  also 
because  the  convergence  of  the  algorithm  is  easier  if  all  the 
nodes  are  connected  to  ground  with  a  capacitance.  For  this 
reason,  it  is  useful  to  define  the  areas  of  the  source  and 
drain  diffusions  even  if  the  layout  has  not  been  done  when 
the  designer  is  writing  down  the  input  for  SPICE.  The  use  of 
the  default  values  DEFAD  and  DEFAS  in  the  .OPTIONS  card  to 
avoid  dangerous  floating  nodes  is  recommended. 
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4.3.9  Equations  for  the  Level  3  Model 

This   model  is  developed  to  simulate  short-ohannel   MOSFET 
which  has  a  channel  length  up  to  2um. 

The  basic  model  of  level  3 

Many  of  the  equations  which  have  been  used  in  this  model 
are  empirical.  The  purpose  is  to  improve  the  precision  and  to 
limit  the  complexity  of  the  calculations.  This  approximation 
results  in  the  development  of  more  manageable  basic  equations 
than  used  in  the  level  2  model.  The  current  equation  in  the 
linear  region  is 

(4. 53)  IDS=IDS[VGS-VTH-(1+FB)/2*VDS]*VDS 

where   FB=(GAMMA*FS)/[2*SQRT(PHI-VBS) ] +FN.   The   FS    is  a 
correction   factor   of   short   channel  effect   and   FN   is   a 

correction   factor  of  narrow  channel  effect.  The  details  of 

correction    factors   and   empirical  data   can   be   found  in 
Reference!!  10]  ■ 

The  dependence  of  mobility  on  the  gate  electric  field  is 
simulated  with  a  simpler  equation  than  that  used  for  the 
level  2  model.  Thus, 

(4.54)  Us=UO/[l+THETA*(VGS-VTH)] 

where  THETA  is  the  SPICE  parameter  for  mobility  modulation. 
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Effect  of  channel  length  and  width  on  threshold  voltage 

Based  on  the  hypotheses  mentioned  in  the  level  2  model, 
the  effect  of  the  different  thickness  of  the  depleted  regions 
in  the  cylindrical  junction  region  and  under  the  channel  has 
been  considered  in  the  calculation  of  the  fixed  charge  QB  in 
this  model.  Then,  The  equation  for  the  threshold  voltage  is 

(4.57)  VTH=VFB+PHI-SIGMA*VDS+GAMMA*FS*SQKT(PHI*-VBS) 

+FN(PHI-VBS) 

where  the  parameter  SIGMA  expresses  empirically  the 
effect  of  the  threshold  voltage  on  VDS.  It  is  known  that  the 
shorter  channel  length  is,  the  lower  the  threshold  voltage. 
This  effect  has  a  linear  relationship  between  the  variation 
of  VTH  and  VDS,  and  allow  us  to  introduce  the  term  SIGMA*VDS. 
The  relationship  between  SIGMA  and  Leff  is  also  empirical: 

(4.58)  SIGMA=ETA*(8. 15e-22 ) /( Cox*Lef f "3 ) 

where  ETA  is  a  static  feedback  on  the  threshold,  and  its 
typical  value  is  1. 

Channel  length  modulation  in  the  saturation  region 

The  equation  of  the  level  3  model  simulates  the   channel 
length  modulation  in  saturation  as 

(4.59)  Leff-L'=SQRT[(EP*XD"2/2)+K*XD-2*(VDS-VD(sat)] 

-EP*XD"2/2 
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where  EP=ID( sat) /[GD(sat ) *Lef f ]  which  is  the  lateral  field  at 
channel  pinch-off  point;  ID(sat)  and  GD(sat)  are  the  drain 
current  and  the  drain  conductance  at  saturation  respectively; 
K  is  an  empirical  fitting  parameter,  which  is  used  to  prevent 
a  negative  channel  length,  with  a  typical  value  of  1. 

Effect  of  channel  length  on  velocity  saturation  of  hot 
electrons 

The  saturation  of  hot  electron  velocity  lowers  the 
conduction  current  in  the  linear  region.  This  model  includes 
a  decrease  in  the  effective  mobility  with  the  average 
electrical  field  between  the  source  and  the  drain[10], 

(4.60)   Ueff=Us/[l+(Us*VDS)/(VMAX*Leff ) 

where  Us  is  the  surface  mobility.  If  the  parameter  VMAX  is 
not  input  by  the  user,  Ueff  is  set  to  Us  and  then  this  hot 
electron  effect  is  not  modeled. 

4.3.10  The  Level  4  Model  in  SPICE 

It  is  also  called  Berkeley  Short-channel  IGFET  Model 
(BSIM)  which  is  built  upon  the  Compact  Short-channel  IGFET 
Model  (CSIM).  Since  a  fully  physics-oriented  modeling 
approach  usually  makes  parameter  extraction  particularly 
difficult,  a  semi-empirical  approach  is  used  to  cope  with  the 
rapid   advances   in   technology  and   to   extract  parameters 
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automatically.  Proc2Mod  is  the  program  which  converts  a 
process  runsheet  file  into  a  sequence  of  .MODEL  cards 
suitable  for  a  SPICE  input  deck.  This  model  has  a  large 
number  of  empirical  parameters,  the  extra  care  has  to  be 
taken  to  define  the  values  of  parameters. 

4.3.11  Table  Lookup  Model 

This  is  not  a  SPICE  model,  but  it  has  been  used 
extensively.  It  represents  the  device  characteristics  in 
terms  of  tables  of  data  points  instead  of  using  analytical 
equations.  These  data  points  can  be  obtained  from 
measurements  on  the  test  devices  or  from  device  analysis 
programs  if  the  fabricated  device  is  not  available.  With  the 
use  of  table  models,  the  model  evaluation  time  can  be  reduced 
considerably  while  a  good  accuracy  is  still  maintained.  Also, 
as  the  technology  progresses,  the  present  models  may  not  be 
adequate  to  characterize  new  devices. 

For  example,  the  drain  current  is  a  function  of  five 
variables:  VGS,  VDS,  VBS,  Leff,  and  Weff.  Modeling  the  drain 
current  as  a  function  of  all  these  variables  will  lead  to  the 
use  of  a  five-dimensional  table,  which  is  impractical. 
Therefore,  knowledge  of  the  device  behavior  is  used  to  reduce 
this  five-dimensional  table.  The  threshold  voltage  is 
principally  a  function  of  VBS.  The  second  order  effects  which 
affect   the  threshold  voltage  are  the  short   channel   effect. 
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the  narrow  channel  effect,  and  the  drain  induced  barrier 
lowering.  Thus,  the  threshold  voltage  is  modeled  by  three 
tables : 

(4.61)  VTH=VT1(VBS)+VT2(VDS,  Wef f )+VT3( VDS,  Leff) 

Then,  a  new  variable  VA  =  VTH  -  VGS  can  be  made.  The  drain 
current  can  then  be  modeled  by  one  three-dimensional  table: 

(4.62)  ID=ID3(VA,  VDS,  VBS) 

Therefore,  the  measured  or  simulated  data  can  be  easily 
transformed  into  the  needed  tables  and  table  dimensions  can 
be  chosen  based  on  sensitivities  of  device  characteristics. 
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4.4   Charge-Control  Model  of  Dynamic  Operation  for   SPICE   in 
Transient  Analysis 

In  the  previous  section,  the  MOSFET  has  been  treated  as 
if  all  terminal  voltage  were  constant.  However,  the  device  is 
usually  used  with  varying  terminal  voltages  in  a  circuit. 
The  varying  terminal  voltages  cause  the  transistor  to  vary 
the  charges  which  must  be  supplied  from  outside  other  than 
transistor  itself.  The  charges  are  supplied  through  the 
device  terminal  by  the  extra  current  which  is  not  predicted 
by  dc  theory.  In  this  section  the  charges  and  terminal 
currents  will  be  evaluated. 

4.4.1  Quasi-static   operation 

With  the  charge-control  model  to  analyze  the  dynamic 
operation,  it  is  assumed  that  the  variation  of  the  terminal 
voltages  is  sufficiently  slow,  so  that  the  device  operates 
quasi-statically.  Therefore,  the  charges  per  unit  area  at 
any  time  are  assumed  identical  to  those  that  would  be  found 
if  dc  voltage  were  used.  The  total  charges  of  each  region 
can  be  given  by 

(4.63)  qG{t>=fG(vI)(t)JvG(t),vB(t>,vs(t)) 

(4.64)  qB(t)=fB(vD(t), vQ(t), vB(t), vs(t)) 

(4.65)  qI(t)=fI(vD(t),vG(t),vB(t),vs(t)) 
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Where  qg  and  qg  are  the  transient  charges  stored  in  the  gate 
and  bulk  region,  respectively,  qj  is  due  to  electrons  in  the 
inversion  (channel)  layer.  These  electrons,  which 
continuously  enter  through  the  source  and  eventually  leave 
the  drain,  are  not  really  stored  in  the  device.  Hence,  q-r  is 
simply  the  total  charge  of  the  electrons  that  happen  to  be  in 
the  inversion  layer  at  any  given  instant. 

Assuming  there  is  no  leakage  in  the  gate,  all  the  gate 
current  attributed  to  the  displacement  current  is  associated 
with  a  changing  gate  charge: 

(4.66)  iG(t)=dqQ/dt 

Similarly,  all  bulk  current  ascribed  to  depleting  or 
recovering  acceptor  atoms  is  associated  with  carriers 
changing   the  charge  in  the  bulk: 

(4.67)  iB(t)=dqB/dt 

Finally,  the  sum  of  the  drain  current  and  source  currents 
represents  the  total  charge  changing  in  the  inversion   layer: 

(4.68)  iD(t)+is(t)=dqj/dt 

Where  ij>(t)  is  a  measure  of  the  number  of  electrons  leaving 
the  device  through  the  drain  per  unit  time  and  -ic(t)  is  a 
measure   of   the  number  of   electrons   entering   the  device 
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through  the  source  per  unit  time.  When  vQ  is  increasing,  the 
number  of  electrons  in  the  channel  must  increase.  Therefore, 
the  rate  of  entering  electrons  from  the  source  must  be 
temporarily  greater  than  the  rate  of  removing  electrons  from 
the  drain.  Thus, 

(4.69)  iD(t)=iT(t)+iDD(t) 

(4.70)  is(t)=-iT(t)+iSD(t) 

Where  iT(t)  is  the  "dc-like"  current,  and  iDD(t)  as  well  as 
iSD(t)  are  responsible  for  a  changing  charge  qj  in  the 
channel.  From  the  last  three  equations, 

(4.71)  igpftj+ipptt^dqj/dt 

It  also  can  be  expressed  as 

(4.72)  dqg/dt+dqpA^dqj/dt 

where  a  changing  charges  dqg  in  time  dt  cause  iSD(t)  and  dqD 
in  time  dt  cause  iDD(t).  To  simplify  the  discussion  in  quasi- 
static  operation,  the  charge  changing  which  causes  current  is 
considered  here.  The  components  of  these  currents  will  hold 
for  Kirchhoff's  law  as  well.  Therefore,  the  quasi-static 
approximation  is  built  to  evaluate  the  strong,  moderate,  and 
weak  inversions,  separately,  in  stead  of  a  general  charge 
sheet  model. 
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4.4.2   Non -Quasi -Static  Operation 

It  is  intuitively  clear  that  the  assumption  of  quasi- 
static  operation  will  be  no  longer  valid  if  the  terminal 
voltages  vary  too  fast. 

The  non-quasi-static  model  is  based  on  an  approximate 
solution  to  the  nonlinear  current-continuity  equation  in  the 
channel[8].  The  non-quasi -static  operation  of  the  MOSFET  is 
a  very  difficult  mathematical  analysis.  In  SPICE3A7  the 
analytic  non-quasi-static  long  channel  MOSFET  for  the 
transient  analysis  has  been  implemented.  The  results  of  this 
model  is  suggested  to  check  accuracy  with  the  PISCES  2-D 
device  simulation. 
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5.0  AN  INTERFACE  PROGRAM  BETWEEN  PISCES  AND  SPICE 

The  methodology  to  link  PISCES  with  SPICE  is  to  an 
interface  program  which  can  accept  the  output  data  of  PISCES 
and  then  calculate  it  to  generate  the  model  parameters  for 
SPICE  simulation.  The  flow  of  the  interface  program  is  shown 
in  Fig.  35.  It  consists  of  three  major  parts:  preprocessing, 
parameter  calculation,  and  post  processing.  The  preprocessing 
and  postprocessing  deal  with  the  output  data  of  PISCES  and 
the  preparation  of  model  parameter  for  SPICE  simulation, 
respectively.  The  parameter  calculation  is  partitioned  into 
three  parts:  calculation  of  linear  region  parameters,  the 
calculation  of  saturation  region  parameters,  and  the 
calculation  of  parasitic  resistance  and  capacitance 
parameters. 

5. 1  Parameter  generation 

With  the  doping  profile,  PISCES  can  provide  carrier 
distribution,  impurity  distribution,  internal  electric  field, 
electrostatic  potential,  mobility,  and  internal  current 
distribution.  Since  the  PISCES  has  the  command  card  called 
SOLVE  which  can  instruct  PISCES  to  generate  a  solution  for 
one  or  more  specified  bias  points,  the  parameters  extracted 
from  subthreshold,  linear,  or  saturation  region  can  be 
handled   easily  with   the  bias   statement.    Six   biasing 
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FIG.  35.  ARCHITECTURE  OF  INTERFACE  PROGRAMED. 
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conditions  shown  in  Fig.  36  are  suggested  for  PISCES  to  run 
and  obtain  the  results  which  can  be  used  for  the  input  of 
interface  program.  In  condition  1  of  PISCES  program,  the 
gate,  drain,  source  and  substrate  are  at  aero  potential, 
Hence,  the  bias-independent  parameters  such  as  LD  and  XJ,  and 
parameters  under  equilibrium  condition,  such  as  PB  can  be 
obtained  by  using  impurity  distribution.  In  condition  2, 
the  subthreshold  parameter  NFS,  Surface-fast  state  density, 
can  be  determined  by  the  following  equation 

(5.1)  NFS=[(   VGS/   log(IDS))*q/kT-l]*Cox/q 

Therefore,  NFS  can  be  measured  from  the  slope  of  the  IDS  vs. 
VGS  curve  in  the  subthreshold  region.  In  condition  3,  the 
carrier  concentration  can  be  used  to  determine  zero-bias 
threshold  voltage.  Because  the  inversion  charge  can  be  found 
by  integrating  the  mobile  carrier  density  over  the  channel 
region[17]  and  the  relationship  among  the  channel  inversion 
charge  QI,  VG,  and  VTO  is  [14] 

(5.2)  QI  =  -Cox'  *  (VG  -  VTO) 

where  Cox'  =  Cox  *  W  *  L.  Therefore,  when  QI  is  obtained,  the 
VTO  is  determined  from  the  above  equation 

(5.3)  VTO  =  QI  /  Cox'  +  VG 

In   condition   4,   a  large  bias  is   applied   to   obtain 
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field-dependent  mobility,  so  that  the  mobility  parameter  can 
be  generated,  such  as  UCRIT,  UEXP, and  THETA.  In  condition  5, 
a  substrate  bias  voltages  is  applied  to  obtain  the  body 
effect  coefficient  (GAMMA)  by  using  carrier  distribution 
under  substrate  bias.  In  condition  6,  a  high  drain  voltage  is 
applied  to  facilitate  the  determination  of  the  channel  length 
modulation  coefficient  (LAMBDA),  as  shown  in  Fig,  37. 
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Fig.  37.  The  determination  of  channel-length  modulation 
coefficient[14]. 


5.2  Calculation  and  Measurement  for  MOSFET  Parameter 

The  methods  of  extracting  the  parameters  of  the  first 
three  existing  levels  of  MOS  models  which  are  implemented  in 
SPICE  for  circuit  simulation  is  very  important.  The  accuracy 
of  the  simulation  depends  upon  the  value  of  input  parameters 
and  the  equations  used  in  the  parameter  extraction  program 
must  be  identical  to  those  used  in  the  circuit  simulator. 
Therefore,  it  is  expected  that  the  I-V  characteristics 
predicted  by  the  equations  of  the  model  must  closely  match 
the  measured  data. 

5.2.1  Parameters  of  the  Level  1  Model 

Equations  of  the  level  1  model  were  given  in  Sec.  4.3.4. 
In  this  level,  a  very  large  and  wide  channel  will  be  used  to 
neglect  possible  differences  between  the  nominal  dimensions 
and  the  real  dimensions. 

Measurements  of  VTO  and  KP  in  linear  region 

In  order  to  account  for  the  condition  that  the 
transistor  is  in  a  linear  region,  it  is  necessary  to  carry 
out  the  measurements  with  a  small  VDS.  Fig.  38  shows  the 
setup  for  the  measurement  in  the  linear  region.  Since  the 
term  VDS  /2  can  be  neglected,  the  current  may  be  expressed  as 

(5.4)   IDS  =  BETA*(VGS  -  VTH)*VDS 
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where  BETA  =  KP  *  (Weff  /  Leff).  The  above  equation  expresses 
a  linear  relationship  between  IDS,  VDS,  and  VDS-VTH,  and 
suggests  the  calculation  of  beta  and  VTO  by  linear 
extrapolation  in  which  the  intercept  is  equal  to  VTO  and 
the  slope  equal  to  BETA.  The  input  parameter  for  SPICE,  VTO, 
is  the  threshold  voltage  VTH  with  zero  voltage  between  source 
and  substrate. 

Finally,  to  make  the  best  choice  of  these  two  parameters 
the  VTO  at  the  lowest  current  and  the  KP  at  the  highest 
voltage  are  chosen. 

Measurements  of  VTO  and  KP  in  the  saturation  region 

This  is  an  alternative  method  to  measure  the  same 
parameters  in  the  saturation  region  by  using  setup  shown  in 
Fig.  39.  The  current  in  the  saturation  region  can  be 
described  in  the  following: 

(5.5)   SQRT(IDS)  =  SQRT(BETA)  *  (VGS  -  VTH) 

The  parameters  can  be  carried  out  by  the  graphical 
calculation  which  is  the  same  as  in  the  previous  section.  It 
is  worth  noting  that  the  calculation  of  the  same  parameters 
with  those  two  methods  does  not  give  equal  results.  In  fact, 
the  simulations  with  the  parameters  calculated  in  the  linear 
region  are  good  only  in  the  linear  region,  while  those  with 
the   parameters  calculated  in  the  saturation  region  are  good 
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Fig-    39.  The  effects  of  bias  VSB  between  the  source  and 
back[14]. 
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in  the  saturation  region. 

Calculation  of  the  body  effect 

Figure  36  shows  the  reverse  bias,  VSB,  applied  between 
the  source  and  the  substrate,  which  has  an  effect  on  the 
theshold  voltage.  The  reverse  bias  voltage  reduces  the  free 
carrier  density  in  the  channel  to  increase  the  threshold 
voltage  of  n-channel  transistor.  The  variation  in  threshold 
voltage  can  be  calculated  as  below: 

(5.6)  VTH=GAMMA*[SQRT(2*PB+VSB)-SQRT(2*PB)] 

where  GAMMA  is  defined  as  body-effect  parameter  and  expressed 
as 

(5.7)  GAMMA=SQRT(2*EPSIL0Nsi*q*NA)/Cox 

5.2.2  Parameters  of  the  Level  2  Model 

The  methods  of  calculating  the  parameters  for  the 
equations  of  the  level  2  model  are  based  on  Section  4.3.5. 

Measurements  in  the  linear  region 

The  parameters  VTO,  KP,  GAMMA,  and  PB  have  the  same 
significance  as  in  the  previous  model  and  are  calculated  in 
the  same  way.  The  parameters  of  Eq. (4-29),  UCRIT,  UEXP,  and 
UTRA,  can  be  calculated  by  a  graphical  method  [9],  The  model 
can  be  rewritten  in  the  following: 
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(5.8)  log(BETA'/BETA)=UEXP*[log(Tox*EPSsi*UCRIT)/K- 

1 og ( VGS-VTH-UTRA* VDS ) ] 

where  BETA'  is  the  value  of  BETA  calculated  with  different 
VGS  and  a  constant  VDS.  Since  UTRA  does  not  have  a  great 
influence  on  the  current  calculation  for  this  model,  it  is 
the  best  choice  to  set  UTRA=0.5.  With  a  suitable  computer,  it 
is  best  to  use  a  linear  regression  method  which  provides  the 
parameters  of  the  linear  equation  y=ax+b,  where 
y=log( BETA' /BETA)  and  x=log(VGS  -  VTH  -  VDS/2 ) .  From  this 
expression  UEXP  and  (JCRIT  are  obtained  as  below: 

(5.9)  UEXP  =  -a 

(5.10)  UCRIT  =  (K/Tox) )*10E-(b/a) 

Calculation  of  physical  parameters 

The  parameters  VTO,  KP,  PB  and  GAMMA  are  of  an 
electrical  type  and  are  related  to  the  parameters  of  a 
physical  type:  mobility  U0,  substrate  doping,  NA,  gate  oxide 
thickness  Tox,  and  number  of  surface  state  NSS.  If  the  value 
of  Tox  is  not  already  known  from  direct  measurements,  it  can 
be  calculated  from  a  measurement  of  the  gate  capacitance,  as 
described  in  Section  4.3.7.  Once  Tox  and  Cox  are  known,  the 
calculation  of  the  mobility  and  the  doping  can  be  obtained  by 
inverting  Eqs .  (4.35)  and  (4.37): 
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(5. 11)  UO  =  KP/Cox 

(5.12)  NA=(Cox"2  *  GAMMA" 2)  /  (2  *  q  *  EPSILONs  *  Ksi) 

Once  NA  is  known,  we  can  recalculate  PB  with  Eq. (4.37). 
If  the  value  found  is  different  from  that  used  to  calculate 
GAMMA,  it  is  better  to  repeat  the  calculation  of  Eq.  (5.7). 

The  threshold  voltage  VTO  can  be  calculated  from 
physical-type  parameters,  using  Eq.  (4.41). 

5.2.3  The  Parameters  of  the  Level  3  Model 

"The  equations  for  the  level  3  model  has  been  given  in 
Sec.  4.3.9.  This  model  is  applicable  to  MOS  transistors  with 
channel  length  of  2um. 

Measurements  in  the  linear  region 

For  this  model,  the  equation  for  simplest  model  in  the 
calculation  is  still  true  for  a  small  VDS.  Once.  VTO,  GAMMA, 
MUO,  NA,  and  PHIp  have  been  found,  we  must  find  the 
parameters  for  the  model  of  the  mobility  variation  with  VGS 
which  has  been  discussed  previously.  The  simplest  method 
consists  of  calculating  BETA'  with  at  least  two  different 
values  for  VGS  and  then  finding  the  values  of  KP  and  THETA 
to  be  inserted  in  the  model  by  solving  a  system  of  two 
equations  and  two  unknowns  as  listed  below: 

(5. 13)      THETA=[(BETA1'/BETA2' ) -1 ]/[VGS2-VT0-(BETAl ' - 
BETA2'  XVGS1-VT0)] 
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(5. 14)      KP=(BETA1'*L/Weff > [ 1+THETA( VGS1-VT0) ] 

He  can  also  obtain  a  better  result  with  an  iterative 
method  which  provides  a  great  precision  by  using  a  least- 
square  curve  fitting  with  a  greater  number  of  measurements 
than  in  the  previous  method. 

Parameters   for   the  short-channel  effect   on   the   threshold 
voltage 

The  model  uses  the  two  parameters  for  the  effect  of  the 
short   channel   on   the  threshold  voltage:  Xj   and   ETA. 

The  value  of  Xj  calculated  in  level  2  model  can  be  used 
only  as  an  initial  value  for  a  trial-and-error  procedure  for 
level  3  model. 

The  parameter  ETA,  proportional  to  the  variation  of  the 
threshold  voltage  with  VDS,  is  easily  calculable  from  the 
threshold  measurements  in  the  saturation  region  and  with  VDS 
constant. 

Parameters  for  the  model  of  the  conductance  in  the  saturation 
region 

Again  a  graphical  method  is  used  to  carry  out  some 
control  simulations  and  correct  the  K  value  by  trial  and 
error. 

The  relationship  between  the  value  of  K  and  the  short 
channel  effect  is  obtained  from  Eq.(4.59).  From  this  equation 
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an   approximate   equation  may  be  obtained,   which  gives   the 
value  of  K  from  the  short  channel  in  the  saturation  region: 

(5.15)       K*XD2=  (Leff-L')2/  VDS 

where  Lef f-L' =Lef f [ l-( IDS/ID( sat) ] . 

The  measurements  of  an  output  characteristics  and 
simulations  with  different  values  for  K  have  shown  that  K  can 
be  obtained  to  represent  the  best  result  in  level  3  model. 
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CONCLUSION 

Simulation  is  an  experiment.  Although  this  report  is  a 
tentative  one,  it  is  worth  integrating  this  environment  for 
both  research  and  teaching  at  a  university  of  this  size  like 
Kansaa  State  University. 

Computer-integrated  manufacturing  is  very  important  in 
the  production  of  VLSI  circuit,  and  SUPREM  and  SPICE  have 
been  widely  used  in  the  IC  industry.  If  these  tools  are 
intelligently  integrated,  a  process  variable,  such  as  ion 
implantation  which  affects  threshold  voltage  and  voltage 
gain,  shown  in  Fig.  40,  can  be  throughout  understood. 

In  this  report,  a  comprehensive  material  of  SUPREM, 
PISCES,  and  SPICE  is  reviewed  and  the  calculation  of 
extracting  MOS  parameter  to  build  an  interface  program  is 
proposed.  We  know  the  model  must  be  verified  by  the 
measurement  taken  on  test  device  which  may  be  fabricated  in  a 
different  foundry.  HSPICE  focuses  on  the  needs  of  the 
integrated  circuit  designer  with  17  MOSFET  models  to  suit  all 
tastes  in  industry.  Sometimes,  the  modification  of  the 
equation  inside  the  simulation  package  is  strongly  suggested 
to   fit  the  measured  data  curve,  such  as  HSPICE   and   PSPICE. 

Therefore,  with  the  knowledge  of  Semiconductor  device, 
programming  language  in  C  and  FORTRAN,  and  numerical 
analysis,  an  integrated  simulation  environment  can  be  built. 
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voltage  gain 

240 


5  10  15  20  25 

implantation  energy  (KeV) 


Fig.  40.  Plot  of  the  voltage  gain  vs.  the  boron 
implantation  energy[l]. 
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APPENDIX     A 
Calcualation  Process  for  PISCES 
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*********    PISCES-II       Device      Modeling       Program       ********* 
(Version    3822-C    sun-68881) 

Date    and    time    =    25-Aug-89    14:55:37 


title  pn  diode  (fig  2a) 

$      Steady-state  -forward  bias  analysis  lw/  contact 

*   resistor) 


t; 

. .  mesh 

4 

.  .  * 

5 

.  .  conta 

6 

. .  symb 

7 

. .  metho 

s 

. .  model 

9 

.  .  * 

10 

- .  sol ve 

11 

.  .  1  og 

12 

. .  sol ve 

13 

. -  solve 

14 

,  ..  sol  ve 

15 

.  .  * 

16 

.  .  end 

i  n-pn . msh 

num=l  resis=le5 

newton  carr=2 

rhsnorm  xnorm  autonr 

temp=300  srh  auger  conmob  tldmob  print 

init  ascii  out-f  =pn2aO--r  .  asc 
out-f  =  IV-r.  log 

vstep=0. 1  nsteps=10  elect=l 
vstep=0.2  nsteps=4  elect=l 
vl=2  ascii  out-f  =pn2a2-r .  asc 


***************** 
pn  diode  (-fig  2a) 

***************** 


Grid  read  -from  pn-msh 
Date  code  =  25-Aug-89  14: 


=h  statistics  : 

Total  grid  points  =  ISO 
Total  no.  o-f  triangles  = 
Obtuse  triangles  =     0 


236 
:  0.0X) 


Newton  method       carriers  - 


Constants  : 

Boltzmanns  k  =  8.61700e-05 

charge        =  1.60200B-19 

permittivity  «  8.85400e-14 

Temperature      =  300. 
Thermal  voltage  "0.025851 


Material  data  : 
Num   Type   Rel  F'ermit 
1      1     1 1 . 80 


Semiconductor  data  : 

type  =   1.0000 

rel  permit  =   11.800 

logKKni  )  =   10.  160 

n-mobility  ■   1000.0 

p-mobility  =   500.00 

vsat  =  1 . 03494e+07 

taupO  =  1 . 00000e-07 

taunO  =  1 . 00000e-07 

eg ap (300)  »   1.0800 

egalpha  =  4.73000e-04 

egbeta  =   636.00 

a-ff  inity  =   4.  1700 

egap  ■   1 . 0800 

cnau  =  2.80000e-31 

cpau  ■  9.90000e-32 

An**  =   110.00 

Ap**  =   30.000 

Nc  =  2.80000e+19 

Nv  =  1.04000e+19 

gcb  =   2.0000 

gvb  =   4.0000 

edb  =  4.40000e-02 

eab  =  4 . 50000e-02 


Contacts  : 
Num     Work  fn.      Vsurfn  Vsur-fp         Resistance 

Capaci  tance 

1      0.000     0.000000e+00    0.000000e+00     1 . 000000e+05 
0. 000000e+00 


u: 


2      0.000     0. 000000e+00    0.000000e+00     0.000000e+00 
. 000000e+00 


Model  flags  : 
SRH  recombination  =  t. 
Cone-  dep.  lifetime  -  f 
Auger  recombination  -  t 
Band-gap  narrowing  -  f 
Cone.  dep.  mobility  =  t 
Field  dep.  mobility  =  t 
Impact  Ionization    =  f 


Mobility  model  parameters 
Gsurf    =   1 . 000000e+00 
B . e 1 ec  t  =   2 . 000000e+00 
B.hole   -   1 . 000000e+00 


Solution  for  bias: 

VI  =   0.0000000e+00       V2  =   0. 0000000e+00 


Initial  solution 


:er 

P1 

51  -error 

n— err or 

p-error 

0 

1 

.  2100e-15 

1, 

.2134e-21 

1. 

. 5428e-19 

RHS 

1 

1 

. 1976e+01 

5 

. 2343e+00 

1 

. 197&e+01 

1 

1 

. 2385e-17 

3, 

, 2185e-06 

1, 

.  M37e-05 

RHS 

2 

6 

. 1309e+00 

2, 

. 8829e+06 

1 

. 7463e+06 

2 

2, 

. 9797e-18 

2, 

. 6373e-06 

1, 

. 3722e-06 

RHS 

3 

1, 

. 9094e+00 

4, 

, 8920e+00 

8 

. 6457e+00 

3 

7, 

, 7556e-19 

1. 

, 1680e-06 

8, 

, 5473e-08 

RHS 

4 

2 

, 1 607e-0 1 

3. 

, 5895e+00 

.j. 

, 1 160e+00 

4 

5 

,  9335e-2<S 

5. 

,8685e-09 

3, 

. 5095e-10 

RHS 

5* 

1, 

, 6560e-03 

9. 

, 9886e-01 

1 , 

. 7B64e+00 

5# 

4, 

7580e-23 

2. 

. 5566e-10 

9, 

. 5634e-12 

RHS 

h 

h. 

, 2307e-05 

9. 

, 9991e-01 

5, 

. 6005e+00 

a 

5 . 

0543e-30 

2. 

1961e-15 

3. 

,8121e-17 

RHS 

7* 

3. 

0450e-10 

8. 

7955e-02 

9, 

. 4912e-05 

7* 

5. 

0953e-30 

2. 

0969e-20 

2, 

, 6292e-lB 

RHS 

B* 

4. 

8217e-14 

1. 

5050e-05 

1. 

, 4257e-09 
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8*     5.036Be-30    3.8861e-21    1.4695e-18   RHS 

9*     4.1777e-14    3.7639e-09    4.1549e-14 

9*     5.0368e-30    3.86426-21    i.4695e-18   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)        (amps/micron)         (amps/micron) 
(amps/micron) 

1  0.0000        1.19554B-31  1.46950e-lB 

1.469508-18 

6. 71969e-22 


0.0000        6.71993e-22  -2.34709e-26 


Electrode  Flux  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

(amps/mi  cron) 

1  -4 . 94367e- 18  0 . 00000e+00 

l.46950e~18 

6.71969e-22 


. 42226e-27  0. 00000e+00 


Absolute  convergence  criterion  met  -for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  for  bias  point  =    34.04 
Total  cpu  time  =    39.00 
Solution  written  to  pn2a0-r.asc 


Solution  for  bias: 

VI  =   1 . 0000000e-01       V2  =   0. 0000000e+00 


Previous  solution  used  as  initial  guess 


ter  psi -error 

0  5.0117e-30 
1'  3.8683e+00 

1  6.4099e-29 

2  9.9459e-02 
2  4.3626e-30 
3*  3.6382e-04 
3*  4.0784e-30 
4*  1 . 3329e-06 
4*  4.1280e-30 
5*  2.2513e-08 
5*  4.0938e-30 


n-error 

p -err or 

3. 

8642e-21 

1. 

0000e-06 

RHS 

3. 

8667e+00 

3. 

, 5394e+00 

4. 

0159e-08 

5 . 

6025e-08 

RHS 

8. 

, 2470e+00 

8. 

. 6174e+00 

1. 

0100e-09 

9, 

, 6420e-ll 

RHS 

6. 

- 2437e-02 

4. 

, 1006e-02 

2. 

0222e-ll 

1. 

. 9965e-13 

RHS 

2. 

3995e-03 

1, 

, 3531e-03 

2. 

5742e-13 

9. 

0828e-16 

RHS 

6. 

2228e-05 

3. 

. 2757e-05 

2. 

62!15e-15 

3. 

,7587e-18 

RHS 

6*     i. 61028-10    1.14078-06    6.6901e-07 

6*     4.1382e-30    2.9127e--17    2.1666e-18   RHS 

Electrode      Voltage  Electron   Current     Hole    Current 
Conduction  Current 

(Volts)  (amps/micron)          (amps/micron) 
( amps/mi cron ) 

1  0.1000  1. 65801 e- 19           8.0165Be-15 
8.01674e~15 

2  0.0000     -7.41B86e-15       -5. 9571  be- lb 
B. 014588-13 

Electrode  Flux  Displacement     Current 

Total  Current 

(caul /micron)  (amps/micron) 

( amps/micron ) 

1  -4.943678-1B  0.00000e+00 

8.016748-15 

8. 014588-13 


■4 .  3 1  845e-26  0 .  00000e+00 


Absolute  convergence  criterion  met  -for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
TDtal  cpu  time  -for  bias  point  =    21.30 
Total  cpu  time  «=    61.14 


Solution  for  bias: 

VI  =   2.0000000e-01       V2  *   0. 0000000e+00 


F'rojection  used  to  -find  initial  guess 

p-error 
6.0842b- 15   RHS 
6.7cb46e+00 
1.5959e-10   RHS 
3 . 0763e-02 
8.39178-13   RHS 
1 . 34000-04 
3.9930(9-16   RHS 
4 . 20868-07 
6.11318-19   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 


er 

psi-error 

i -error 

0 

7. 

2841e-19 

1, 

. 3846e-14 

1 

2 

07568-01 

3, 

. 9941e-01 

1 

4. 

2797e-30 

7, 

, 20218-09 

2 

1, 

1284e-02 

1, 

, 38968-02 

2 

4. 

1 559e-30 

1, 

,5215e-U 

3  » 

4. 

8082e-06 

3, 

, 6296e-05 

3* 

4. 

10298-30 

2, 

, 8475e-14 

4* 

2, 

05438-09 

3. 

0406e~07 

4* 

4. 

10028-30 

4, 

, 49108-17 

(amps/mi  cron) 

1  0 . 2000        8.101 30e~ IS  9 . 48662e- 1 4 
9. 48743e- 14 

2  0.0000     -6. 649928-14       -2.837568-14 

9.4B747e-14 

Electrode  Flux  Displacement     Current 

Total  Current 

(caul /micron)  (amps/micron) 

(amps/micron) 

1  -4.94367e~18  0.00000e+00 

9.48743e-14 

9.487478-14 


-7.0941 8e— 25  0 .  00000e+00 


Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  -for  bias  point  =    15.90 
Total  cpu  time  "    77.08 


Solution  -for  bias: 

VI  -      3. 0000000e-01       V2  =   0. 0000000e+00 


Projection  used  to  -find  initial  guess 

iter     psi -error  n— error       p—srror 

0  B. 820ie-19  3.2788e-14    1.1622e-13   RHS 

1  2.4027«-01  2.40288-01    2.3596e-01 

1  4.3728e-31  1.0156e-08    1.9461e-10   RHS 

2  1.5993e-02  1.91868-02    4. 1460e-02 

2      1.73748-30  3.2754e-ll    1.8115e-12   RHS 

3*     9.8197e-06  7.07758-05    2.6116e-04 

3*     3.7331e-30  8.7661e-14    1.51468-15   RHS 

4*     6.2684e-09  8.2234e-07    1.15568-06 

4*     3.74B0e-30  1.50568-16    8.63608-19   RHS 

Electrode      Voltage  Electron   Current     Hole    Current 
Conduction  Current 

(Volts)  (amps/micron)         (amps/micron) 
( amps/mi  cron) 

1  0 . 3000  3 . B7903e- 16           1 . 94954b- 1 2 
1.94993e-12 

2  0.0000  -6.29432e-13       -1.320508-12 
1.94993e-12 
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Electrode  Flu::  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

(amps/micron) 

1  -4.94367e-18  0.00000e+00 

1.94993e-12 

1.94993e-12 


:>6944e--23  0. 00Q00e+00 


Absolute  convergence  criterion  met  -for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  -for  bias  point  =    16.02 
Total  cpu  time  =    93.14 


Solution  far  bias: 

VI  =   4. 000000 le— 01       V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 

iter     psi -error"      n— error       p-error 

0  1.1009e-18    2. 966 le- 13    3. 225 le- 12   RHS 

1  2.8302e-01    2. 8306e-01    2.7347e-01 

1  3.6992e-3o    1 . 6990e-08    3.4661e-10   RHS 

2  2.3163e-02    2.69128-02    5.8402e-02 

2  6.63738-30  7.7554e-ll  4.34048-12   RHS 

3*  2.16328-05  1 . 5627e-04  5.3143e-04 

3*  4.8127e-30  3. IOSSb-13  6.0682e-15   RHS 

4*  2.2493e-08  2.5906e-06  3.4023e-06 

4*  2.63588-30  9.9174e-16  2.2827e-17   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/micron) 

1  0.4000        1.85600e-14  6.95404e-ll 
6.955908-11 

2  0.0000     -B.  51892e--12       -6.10401e-ll 
6.955908-11 

Electrode  Flu::  Displacement     Current 

Total  Current 

(coul /mi cron)  (amps/micron) 

(amps/micron) 

1  -4 . 94366e- 18  0. 00000e+00 

6.95590e-ll 
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2  -8.53301e— 22         0.00000e+00 

6. 95590e-ll 

Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  for  bias  point  =    15.98 
Total  cpu  time  =   109.18 


Solution  for  bias; 

VI  =   5. 000000 le-01       V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 

iter     psi -error      n -error       p-error 

0  1.4775e-18    3.3138e-12    i.4817e-10   RHS 

1  3.0428e-01    3.0516e-01    2. 741 le-01 

1  5.3512e-30    3.14376-08    1 . 0043e-09   RHS 

2  2.4519e-02    3.5413e-02    h. 1714e-02 

2  8.6700e-31  2.1464e-10  1.8389e-ll   RHS 

3*  2.3502e-05  2.6529e-04  5.9717e-04 

3*  5.253Se-30  9. 653 le- 13  1.0223e-13   RHS 

4*  2.9B70e-07  3.6458e-06  2.2941e-06 

4*  1. 30920-30  3.61916-15  4.9477e-16   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

( amps/mi cran) 

1  0.4997        B.77881e-13  2.87792»-09 
2 . 87880e -09 

2  0 . 0000     -1.91 7&7e- 1 0       -2 . 68703e-09 
2.87880e-09 

Electrode  Flux  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

(amps/mi  cr on ) 

1  -4.94363e-18  0.00000e+00 

2.87880e-09 

2.87880e-09 


-  3 .  6  4  /  0  9  e  -  2  0  0 .  0  0  0  0  0  e + 0  0 


Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  for  bias  point  =    16.20 


Total  cpu  ti me  =   125. 41 


Solution  -for  bias: 

V 1  ~   6 .  0  0  0  0  0  0 1  e  -  0 1       V  2  -   0 .  0  0  0  0  0  0  Oe + 0  0 


Projection  used  to  find  initial  guess 


er 

psi —error 

n -err or 

p— error 

0 

1. 

, 3523e-17 

6. 

5537e-ll 

5. 

6043e-09 

RHS 

1 

4, 

, 2455e-01 

4. 

1395e-01 

6. 

, 8982e-01 

1 

3, 

, 5192e-30 

1. 

1469e-07 

2. 

. 8199e-08 

RHS 

2 

7, 

, B559B-02 

1. 

5505e-0i 

4 . 

, 5004e-01 

2 

6, 

, 2293e-30 

2. 

2861e-09 

8. 

7123e-10 

RHS 

3  * 

8, 

, 6780e-03 

2. 

9395e-03 

7. 

, 3773e-03 

3* 

Si 

,6791e-30 

3. 

3314a- 10 

1. 

4653e-10 

RHS 

4 

1. 

, 19366-03 

4. 

0592e-04 

1. 

, 0500e-03 

4 

9. 

, 5052e-31 

1. 

0933e-13 

D< 

7105e-14 

RHS 

5* 

4. 

, 0932e-07 

1. 

1519e-07 

3. 

, 2254e-07 

5* 

4, 

, 1612e-30 

7. 

<b794e-17 

. 7435e-17 

RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/ micron)  (amps/micron) 

(amps/rni  cr on) 

1  0.5934        3. 2592 le- 11  6.63809e-08 
6. 64135e-08 

2  0.0000     -2.53824e-09       -6. 38752e-08 
6.64135e-08 

Electrode  Flux  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

( amps /micron) 

1  -4.942B7e-18  0.00000e+00 

6. 64135e-08 

6. 64135e-08 


-7. B0089e-19  0. 00000e+00 


Absolute  convergence  criterion  met  -for  Poisson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  for  bias  point  «    20.28 
Total  cpu  time  =   145.76 
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Solution    for    bias: 

VI  =   7. 000000 le-01       V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 

iter     psi-error      n— error       p-error 

0  1.1577e-16    2.9<b71e-10    5.5512e-08   RHS 

1  7.8065&-01    7.9635e-0i    B.5B14e-01 

1  1.4157e-29    1.2302e-G6    2. 6B40e-07   RHS 

2  4.4696e-01    4.3923e-01    5.8549e-01 

2  2.7533e-30    1 . 33S9e-07    3.4278e-08   RHS 

3  3.8393e-02    3. 2948e-02    1.3314e-02 

3  7.1658e-30  6.2345e-10  4.2224e-10   RHS 

4*  1.8952e-04  9. 6796e-05  5.638ie-05 

4»  3.1167e-30  1.1933e--12  5.0417e-12   RHS 

5*  1.5970e-06  1 . 2470e-06  6.5157e-07 

5*  6.4227e-30  3.0605e-14  1.2541e-14   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/micron) 

1  0.6616        4.50574e-10  3.83823e-07 
3.B4273e-07 

2  0.0000     -1.01286e-08       -3.74145e-07 
3.84273e-0  7 

Electrode  Flux  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

( amps/mi  cr on) 

1  -4.93904e-lB  0.00000e+00 

3.84273e-07 

3.842  73e-07 


3 . 4004 1 e- 1 8  0 . 00000e+00 


Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  for  bias  point  =    20. 2B 
Total  cpu  time  =   166.08 


Solution  for  bias: 

VI  =   8. 000000 le-01       V2  =   0. 0000000e+00 


1  2 1 


Projection  used  to  find  initial  guess 

iter     psi -error      n— error  p-error 

0  5.9720e-17    1 . 9888e-09    6. Q837e-08   RHS 

1  8.0458e-01    8.0437e-01    5.9B24e-01 

1  3.1514e-30    2.7519e-06    1 . 37088-07   RHS 

2  3.6788e 01    6.  7998e-01    1.2635e-01 

2      2.7670e-30  3.8219e-07  1 . 3B34e-07   RHS 

9.0363e-03  3.9647e-02  1 . 24B5e-03 

5.7398e-30  7.3493e-10  1.5573e-10   RHS 

4*     9.2535e-06  4. 10238-05  4.0086e-06 

4*     6.B047S-30  7.4215e-13  2. 5812e-13   RHS 

5*     3.8745e--09  7.5614e-08  7.8460e-09 

5*     4.6911e-30  2.4062e-15  9. 1422e-16   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/micron ) 

1  0.7033        2.23499e-09  9.65054e-07 
9.67289e-07 

2  0.0000     -2.32936e-08       -9. 43996B-07 
9.67289e-07 

Electrode  Flux  Displacement     Current 

Total  Current 

Ccoul /micron)  (amps/micron) 

( amps/mi  cron ) 

1  -4 . 93203e- 1 S  0 . 00000e+00 

9.67289e-07 

9.67289e-07 


-6 . 35497e- 1 B  0 . 00000e+00 


Absolute  convergence  criterion  met  for  Poi sson 
Absolute  convergence  criterion  met  for  continuity 
Total  cpu  time  -for  bias  point  =    20.16 
Total  cpu  time  =   186.30 


Solution  for  bias: 

VI  =   9. 000000 le-01        V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 
iter     psi -error      n -err or       p-error 
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0 

3. 

,6816e-17 

3. 

,  0777e- 

■09 

1 

5. 

2005e-01 

5. 

1904e- 

-01 

1 

2, 

. 9260e-30 

1. 

8243e- 

-06 

2 

9, 

8961e-02 

1. 

.  3929e- 

-01 

2 

4, 

9049e-30 

3 . 

963Be- 

-08 

'  * 

B. 

9552e-04 

2. 

6231e- 

-03 

3* 

!■  . 

5521e-30 

2, 

1093e- 

-09 

4# 

1. 

727ie-05 

3. 

1759e- 

-05 

4* 

2. 

7282e-30 

6. 

3730e- 

-11 

5* 

3. 

2001e-07 

9. 

0229b- 

-07 

5* 

7, 

9250e-31 

1. 

7693e- 

-12 

B.3561e~08   RHS 
2.B630e-01 
1.6149e-07   RHS 
3. 1891e-02 
2.2731e-08   RHS 
7.9452e-05 
7.417Be-10   RHS 
1.7214e-06 
2.3039e-ll   RHS 
5.3047e-08 
6.5175e- 13   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/mi  cron ) 

1  0.7303        6.28149e-09  1 . 69093e-06 
1.69721e-06 

2  0.0000     ~4. 10851 e-08       -1.65612e-06 
1.69721e-06 

Electrode  Flux  Displacement     Current 

Total  Current 

(caul /micron)  (amps/micron) 

(amps/mi  cron ) 
1  -4.92324e-18  0.00000e+00  1.69721e-06 

2  -8.79339e-18  0.00000e+00  - 

1.69721e-06 

Absolute  convergence  criterion  met  -for  Poisson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  -for  bias  point  =    18.74 
Total  cpu  time  =   205.08 


Solution  -for  bias: 

VI  =   1 . 0000000e+00       V2  =   0. 0000000e+00 


Projection  used  to  -find  initial  guess 

iter  psi— error  n— error  p-error 

0  2.4813e-17  2. 9204e-09  7.0269e-08   RHS 

1  2.9356e-01  2.9339e-01  1 . 4489e-01 

1  9.1062e-30  7.6141e-07  1 . 0034e-07   RHS 

2  2.6091e-02  3. 0830e-02  7.0677e-03 

2  5.2757e-30  5.0476a-09  2.1278e-09   RHS 
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3*  4.9525e-05  1.3049e-04  7.3957e-06 

3*  8.B796e-31  4. 276 ie- 11  1.8058e-ll   RHS 

4*  2.6784e-07  3.8261e-07  3.6610e-0S 

4*  1.6224e-30  3.2709e-13  1.4107B-13   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/mi cron ) 

1  0.7496        1.30938e--08  2.49122e-06 
2.50431e-06 

2  0.0000     -6. 29148e-08       -2.44140e-06 
2.50431e-06 

Electrode?  Flux  Displacement     Current 

Total  Current 

(coul /micron )  (amps/micron) 

(amps/micron) 

1  -4.91 354e- 18  0 . 00000e+00 

2. 50431e-06 

2. 50431e-06 


- 1 . 07556e- 1 /  0 . 00000e+00 


Absolute  convergence  criterion  met  -for  Poi sson 
Absolute  convergence  criterion  met  tor  continuity 
Total  cpu  time  for  bias  point  =    16.10 
Total  cpu  time  =   221.24 


Solution  tor  bias: 

VI  =   1.2000000(8+00       V2  =   0. 0000000e+00 


Projection  used  to  -find  initial  guess 


p -err or 
1.8208e-07   RHS 
2.0778e-01 
3.7981e-07   RHS 
1 -5431e-02 
1 . 4206e-08   RHS 
4.5482e-05 
2.5578e-10   RHS 
5. 8064e-07 
4.2371e-12   RHS 


iter 

psi -error 

r 

i-error 

0 

7. 

. 8735e-17 

9. 

2779e-09 

1 

4. 

. 2213e-01 

4, 

2172e-01 

1 

2, 

, 2135e-29 

2 . 

4853e-06 

2 

6, 

, 0212e-02 

7, 

6151e-02 

2 

3 

. 5248e-30 

2. 

S580e-08 

3* 

2, 

, 2689e-04 

b. 

8862e-04 

3* 

2. 

2829e-30 

4. 

6016e-10 

4* 

3. 

1132e-06 

5. 

2624e-06 

4* 

3. 

8809e-30 

7 . 

5135e-12 

EI ec trod 

e 

Vol tage 

Electroi 
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Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/micron) 

1  0.7764         3.60073e-08  4.20039e-06 

4. 23640B-06 
2         0.0000     -1.17835e-07       -4. 118566-06      -4.23640e-06 

Electrode  Flux  Displacement     Current 

Total  Current 

(caul /mi cron )         (amps/micran)  (amps/micran; 

1  -4 . 3927 1 e- 1 8  0 . 00000e+00 

4.23640e~06 
2  - 1 . 37459s- 17  0 .  00000e+00  -4 . 23640e-06 

Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  for  bias  point  -    16.14 
Total  cpu  time  =   237.50 


Solution  for  bias: 

VI  =   1 . 4000000e+00       V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 


er 

psi -error 

n-error 

p -err or 

0 

7. 

, 0956e-17 

1. 

1299e-08 

1. 

, 8166e-07 

RHS 

1 

3, 

,0316e-01 

3. 

, 0408e-01 

1, 

, 466Se-01 

1 

2, 

, 4420e-30 

1. 

7799e~06 

2, 

, 3825e-07 

RHS 

2 

2, 

, 897Be-02 

3. 

. 1728e-02 

7, 

,  1459e-03 

2 

6. 

. 8119e-30 

1. 

1 989e-08 

5. 

1818e-09 

RHS 

3* 

4. 

. 7B17e-05 

1, 

1999e-04 

1. 

, 2226e-05 

3* 

2. 

3695e-30 

1. 

1007e-10 

4. 

6882e-ll 

RHS 

4* 

3. 

, 6759e-07 

3, 

8623e-07 

7, 

, 9503e-08 

4* 

2, 

4214e-30 

9. 

1949e-13 

4. 

0002s- 13 

RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/micron) 

1  0.7950        7.21161e-08  5.97770e-06 
6. 04982e-06 

2  0.0000     ~1.87292e-07       -5.B6252e-06 
6. 04982e-06 

Electrode  Flux  Displacement     Current 
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Total  Current 

(coul /micron)  (amps/micron) 

(amps/mi  cron) 

1  -4 . B709 1 e- 1 8  0. 00000e+00 
6.04982e-06 

2  -1.5991 0e- 17  0.00000e+00 
6. 04982e-06 

Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  for  bias  point  =    16.14 
Total  cpu  time  =   253.68 


Solution  for  bias: 

VI  =   1 .  6000000e+00       V2  =   0. 0000000e+00 


Projection  used  to  find  initial  guess 

iter     psi— error      n—error  p-error 

0  5.0683e-17    9.9209e-09    1 . 2544e-07   RHS 

1  1.6620e-01    1.6811e-01    8.0620e-02 

1  7.6434e-30    4.3449e-07    1.1405e-07   RHS 

2  7.8062e-03    8.3218e-03    1 . 5496e-03 

2      6.2797e-30    6.5417e-10    2.4498e-10   RHS 

3*     2.3420e-06    7.2072e-06    8.3775e-07 

3*     5.0102e-30    1.2884e-12    4.4248e-13   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)      (amps/micron)       (amps/micron)      (amps/micron) 

1  0.8094        1.21876e-07  7.78460e-06 
7.90647e-06 

2  0.0000     -2.71059e-07       -7.63541e-06 
7.90647e-06 

Electrode  FIuk  Displacement     Current 

Total  Current 

(coul /mi cron)         (amps/micron)  (amps/micron) 

1  -4.84860e-18  0. OOOOOe+00 
7.90647e-06 

2  -1.78010e-17  0.00000e+00  - 
7. 90647e-06 

Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  for  continuity 
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Total  cpu  time  -for  bias  point  ~    13.46 
Total  cpu  time  =   267.18 


Solution  for  bias: 

VI  =   1 . 8000000e+00       V2  =   0. 0000000e+00 


Projection  used  to  -find  initial  guess 

iter     psi -error      n— error       p-error 

0  3.70356-17    9.24218-09    9.6201e-08   RHS 

1  1 . 0285e-0 1    1 .  OS  1 1 B-0 1    5 . 0297e-02 

1  3. 19366-30    2.2742e-07    4.9843e-08   RHS 

2  3.0035e-03    3.0520e-03    5.9616e-04 

2      7.1 498e-30    1 . 0722e- 10    5 . 37 1 1 e- 1 1   RHS 

3*     3.4331e-07    8.4649e-07    1.1504e-07 

3*     7.1456e-30    8.2466e--14    4.36976-14   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)        (amps/micron)         (amps/micron) 
(amps/micron) 

1  0.8210        1.85455e-07  9.604126-06 
9.789586-06 

2  0.0000     -3.689826-07       -9. 42059e-06 
9. 78958e-06 

Electrode  Flux  Displacement     Current 

Total  Current 

(coul /mi cron)  (amps/micron) 

( amps/mi  cron) 

1  -4 . 82599e- 18  0 . 00000e+00 

9.789586-06 

9.789586-06 


- 1 . 93326e- 1 7  0 . 00000e+00 


Absolute  convergence  criterion  met  for  Poisson 
Absolute  convergence  criterion  met  tor  continuity 
Total  cpu  time  -for  bias  point  =    13.46 
Total  cpu  time  =   280.70 


Solution  -for  bias: 

VI  =   2. 0000000e+00       V2  =   0. 0000000e+00 
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Projection  used  to  find  initial  guess 

iter     psi -error      n-error       p-error 

0  2.8067e-17    B.8739e-09    7.2211e-08   RHS 

1  6.9579e-02    7. 1995e-02    3.4160e-02 

1  4.5438e-30    i.2321e-07    2.0-340e-08   RHS 

2  1 . 3432e-03    i . 354Be-03    2.6637e-04 

2      3.1803e-30    2.9961e-li    1.3607e-ll   RHS 

3*     6 . 8924e-08    1 . 5020e-07    2.71 68e-08 

3*     2.6803e-30    i.0378e-14    5.0536e-15   RHS 

Electrode      Voltage     Electron   Current     Hole    Current 
Conduction  Current 

(Volts)         (amps/micron)  (amps/micron) 

(amps/mi cron  > 

1  0.8310        2.62872e-07  1.14275e-05 
1 . 1 6903e-05 

2  0.0000     -4.80902e-07       -1.12094e-05 
1 . 16903e-05 

Electrode  Flux  Displacement     Current 

Total  Current 

(coul /micron)  (amps/micron) 

(amps/mi  cron) 

1  -4.80318e-18  0.00000e+00 

1. 16903e-05 

1 . 16903e-05 


-2 . 06744e- 17  0 . 00000e+00 


Absolute  convergence  criterion  met  tor  F'oi  sson 
Absolute  convergence  criterion  met  -for  continuity 
Total  cpu  time  -for  bias  point  =    13.48 
Total  cpu  time  ■   294.30 
Solution  written  to  pn2a2-r.asc 
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APPENDIX     B 


MOSFET  Model  Parameters  for  SPICE3A7[6] . 
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parameter 


example 


1 

LEVEL 

2 

VTO 

3 

KP 

4 

GAMMA 

5 

PHI 

6 

LAMBDA 

7 

RD 

8 

RS 

9 

CBD 

10 

CBS 

11 

IS 

12 

PS 

13 

CGSO 

14 

CGDO 

15 

CGBO 

16 

RSH 

17 

a 

IS 

Ml 

19 

CJSW 

20 

MJSW 

TON 

NSUB 

NSS 

NFS 
TPC 


27  XJ 

28  LD 

29  UO 

30  UCRIT 

31  UEXP 

32  LTTRA 

33  VMAX 

34  NEFF 


model  index 

zero-bias  threshold  voltage 

transconduciance  parameter 

bulk  threshold  parameter 

surface  potential 

channel-length  modulation 

CMOS!  and  MOS2  only) 

drain  ohmic  resistance 

source  ohmic  resistance 

2ero-bias  B-D  junction  capacitance 

zero-bias  B-S  junction  capacitance 

bulk  junction  saturation  current 

bulk  junction  potential 

gate-source  overlap  capacitance 

per  meter  channel  width 

gate-drain  overlap  capacitance 

per  meter  channel  width 

gate-bulk  overlap  capacitance 

per  meter  channel  length 

drain  and  source  diffusion 

sbeei  resistance 

zero-bias  bulk  junction  bottom  cap. 

per  sq-meter  of  junction  area 

bulk  junction  bottom  grading  coef. 

zero-bias  bulk  junction  sidewall  cap. 

per  meter  of  junciion  perimeter 

bulk  junction  sidewall  grading  coef. 

bulk  junction  saturation  current 
per  sq-meter  of  junction  area 
oxide  thickness 
substrate  doping 
surface  stale  density 
fasi  surface  state  density 
type  of  gate  material: 

+1  opp.  to  substrate 

-1  same  as  substrate 
0   A I  gate 
metallurgical  junction  depth 
lateral  diffusion 
surface  mobility 
critical  field  for  mobility 
degradation  (MOS2  only) 
critical  field  exponent  in 
mobility  degradation  (MOS2  only) 
transverse  field  coef  (mobility) 
(deleted  for  MOS2) 
maximum  drift  velocity  of  carriers 
total  channel  charge  (fixed  and 
mobile)  coefficient  (MOS2  only)- 


- 

1 

V 

0.0 

1.0 

A/V"2 

2.0E-5 

3.1E-5 

V-0.5 

0.0 

0.37 

V 

0.6 

0.65 

1/V 

0.0 

0.02 

Ohm 

0.0 

1.0 

Ohm 

0.0 

1.0 

F 

0.0 

20FF 

F 

0.0 

20FF 

A 

1.0E-14 

1.0E-15 

V 

0.8 

0.S7 

F/a 

0.0 

4.0E-11 

F/bi 

0.0 

4.0E-11 

F/m 

0.0 

2.0E-10 

Ohm/sq. 

0.0 

10.0 

F/m"2 

0.0 

2.0E-4 

- 

OJ 

0.5 

F.  a 

0.0 

1.0E-9 

~ 

O.SOUevell) 
0.330evel2.3) 

A/m"2 

1.0E-8 

meler 

1.0E-7 

1.0E-7 

17cm"3 

0.0 

4.0E15 

lA:m"2 

0  0 

1.0E10 

l/cm"2 

0.0 
1.0 

1.0E10 

meter 

0.0 

1U 

meter 

0.0 

0.8U 

cm«"2/V-s 

600 

700 

V/cm 

1.0E4 

1.0E4 

- 

0.0 

0.1 

- 

0.0 

0.3 

Bfi/l 

0.0 

5.0E4 

1  30 


35 

KF 

36 

AF 

37 

FC 

39 

THETA 

40 

ETA 

41 

KAPPA 

flicker  noise  coefficient 
flicker  noise  eiponenl 
coefficient  for  forward-bias 
depletion  capacitance  formula 
width  effect  on  threshold  voltage 
CMOS2  and  MOS3) 
mobility  modulation  CMOS3  only) 
static  feedback  (MOS3  only) 
saturation  field  factor  (MOS3  only) 


SPICE  BSIM  (level  a)  parameters. 


0.0 
1.0 


0.5 


1.0E-26 
1.2 


0.0 

1.0 

0.0 

0.1 

0.0 

1.0 

0.2 

0.5 

VFB 
PHI 
Kl 
K2 

ETA 
MUZ   ' 
DL 
DW 

uo 

Ul 

X2MZ 

X2E 

X3E 

X2U0 

X3U1 

MUS 

X2MS 

X3MS 

X3U1 

TOX 

TEMP 

VDD 

CGDO 

CGSO 

CGBO 

XPART 

NO 

NB 

ND 

RSH 

JS 

PB 

MJ 

PBSW 

MJSW 

CJ 

CJSW 

WDF 

DELL 


parameter  units 

flat-band  voltage  V 

surface  inversion  potential  V 

body  effect  coefficient  V 
drain/source  depletion  charge  sharing  coefficient 
2ero-bias  drain-induced  barrier  lowering  coefficient 

zero-bias  mobility  cm  /V-s 

shortening  of  channel  M™ 

narrowing  of  channel  Mm 

zero-bias  transverse-field  mobility  degradation  coefficient  V~* 

zero-bias  velocity  saturation  coefficient  fimfV 

sens,  of  mobility  to  substrate  bias  it  v^*  0  cm  /V  -s 

sens,  of  drain-induced  barrier  lowering  effect  to  substrate  bias  V"* 
sens,  of  dnin-induced  barrier  lowering  effect  to  drain  bias  at  V^*  Vw       V* 

sens,  of  transverse  field  mobility  degradation  effect  to  substrate  bias  V"* 

sens,  of  velocity  saturation  effect  to  substrate  bias  /imV 

mobility  at  zero  substrate  bias  and  at  V^s  Vdd  cm2/V!-s 

sens,  of  mobility  to  substrate  bias  at  Vdl*  Vad  cm:/V:-s 

sens,  of  mobility  to  drain  bias  at  V^a  Vdd  cm2/V2-* 

sens,  of  velocity  saturation  effect  on  drain  bias  at  V&=  vta  aud V 

gate  onde  thickness  m-ti 

temperature  at  which  parameters  were  measured  "C 

measurement  bias  range  V 

gate-drain  overlap  capacitance  per  meter  channel  width  F/m 

gate-source  overlap  capacitance  per  meter  channel  width  F/m 

gate-bulk  overlip  capacitance  per  meter  channel  length  F/m 
gate-oxide  capacitance  charge  model  flag 
zero-bias  subthreshold  slope  coefficient 
sens,  of  subthreshold  slope  to  substrate  bias 
sens,  of  subthreshold  slope  to  drain  bias 

drain  and  source  diffusion  sheet  resistance  O/U 

source  drain  junction  current  density  A/m* 

built  in  potential  of  source  drain  junction  V 
Grading  coefficient  of  source  drain  junction 

built  in  potential  of  source.drain  juniion  sidewall  V 

griding  coefficient  of  source  drain  junction  sidewall 

Source  drain  junction  capacitance  per  unit  area  F/m1 

source  drain  junction  sidewall  capacitance  per  unit  length  F/m 

source  drain  junction  default  width  m 

Source  drain  junction  length  reduction  m 
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ABSTRACT 

With  the  increasing  complexity  of  fabrication  processes 
and  decreasing  physical  sizes  of  the  device  structures  in 
integrated  circuits,  integration  of  computer  aided  design 
tools  is  strongly  needed  for  the  automation  of  simulation 
procedures  as  well  as  for  facilitating  the  designers  to 
examine  how  a  fabrication  variable  affects  the  final  device, 
the  circuit  performance,  and  the  product  yield.  This  report 
describes  an  integration  of  process,  device,  and  circuit 
simulators  which  are  SUPREM,  PISCES,  and  SPICE,  respectively. 
Since  the  interface  between  SUPREM  and  PISCES  has  been 
completed,  the  methodology  to  write  an  interface  program 
which  links  a  device  simulator  with  a  circuit  simulator  is 
presented.  Also,  acceptable  MOSFET  models  and  extractions  of 
MOSFET  parameters  for  the  SPICE  simulator  are  reviewed  and 
emphasized. 


